Permit non-conserving initialization for EvoAtmosphere

src/input/photochem_input_after_read.f90

@@ -85,39 +85,17 @@ subroutine interp2atmosfile(dat, var, err)
       enddo
       var%usol_init = 10.0_dp**var%usol_init
 
-    else; block
-      real(dp) :: dz_file
-      real(dp), allocatable :: densities_file(:,:) ! molecules/cm3
-      real(dp), allocatable :: ze_file(:), ze(:)
-
-      dz_file = dat%z_file(2)-dat%z_file(1)
-
-      allocate(densities_file(dat%nq,dat%nzf))
-      allocate(ze_file(dat%nzf+1))
-      allocate(ze(var%nz+1))
-
-      do i = 1,dat%nq
-        densities_file(i,:) = dat%mix_file(i,:)*dat%den_file
-      enddo
-
-      ze_file(1) = dat%z_file(1) - 0.5_dp*dz_file
-      do i = 1,dat%nzf
-        ze_file(i+1) = dat%z_file(i) + 0.5_dp*dz_file
-      enddo
-      ze = var%z(1) - 0.5_dp*var%dz(1)
-      do i = 1,var%nz
-        ze(i+1) = var%z(i) + 0.5_dp*var%dz(i)
-      enddo
+    else
 
-      do i = 1,dat%nq
-        call conserving_rebin(ze_file, densities_file(i,:), ze, var%usol_init(i,:), ierr)
-        if (ierr /= 0) then
-          err = 'subroutine conserving_rebin returned an error'
-          return
-        endif
-      enddo
+      if (dat%conserving_init) then
+        call interp2atmosfile_mixconserving(dat, var, err)
+        if (allocated(err)) return
+      else
+        call interp2atmosfile_mix(dat, var, err)
+        if (allocated(err)) return
+      endif
 
-    endblock; endif
+    endif
 
     if (dat%there_are_particles) then
       do i = 1,dat%npq
@@ -141,6 +119,110 @@ subroutine interp2atmosfile(dat, var, err)
 
   end subroutine
 
+  subroutine interp2atmosfile_mixconserving(dat, var, err)
+    use futils, only: interp, conserving_rebin
+    use photochem_const, only: small_real
+    type(PhotochemData), intent(in) :: dat
+    type(PhotochemVars), intent(inout) :: var
+    character(:), allocatable, intent(out) :: err
+
+    integer :: i, ierr
+    real(dp) :: dz_file
+    real(dp), allocatable :: densities_file(:,:) ! molecules/cm3
+    real(dp), allocatable :: ze_file(:), ze(:)
+
+    dz_file = dat%z_file(2)-dat%z_file(1)
+
+    allocate(densities_file(dat%nq,dat%nzf))
+    allocate(ze_file(dat%nzf+1))
+    allocate(ze(var%nz+1))
+
+    do i = 1,dat%nq
+      densities_file(i,:) = dat%mix_file(i,:)*dat%den_file
+    enddo
+
+    ze_file(1) = dat%z_file(1) - 0.5_dp*dz_file
+    do i = 1,dat%nzf
+      ze_file(i+1) = dat%z_file(i) + 0.5_dp*dz_file
+    enddo
+    ze = var%z(1) - 0.5_dp*var%dz(1)
+    do i = 1,var%nz
+      ze(i+1) = var%z(i) + 0.5_dp*var%dz(i)
+    enddo
+
+    do i = 1,dat%nq
+      call conserving_rebin(ze_file, densities_file(i,:), ze, var%usol_init(i,:), ierr)
+      if (ierr /= 0) then
+        err = 'subroutine conserving_rebin returned an error'
+        return
+      endif
+    enddo
+
+  end subroutine
+
+  subroutine interp2atmosfile_mix(dat, var, err)
+    use futils, only: interp, conserving_rebin
+    use photochem_const, only: small_real
+    type(PhotochemData), intent(in) :: dat
+    type(PhotochemVars), intent(inout) :: var
+    character(:), allocatable, intent(out) :: err
+
+    integer :: i, ierr, nzf1
+    real(dp), allocatable :: z_file(:), den_file(:), density(:)
+    real(dp) :: slope, y_intercept, den_tmp
+
+    z_file = dat%z_file
+    den_file = dat%den_file
+    nzf1 = dat%nzf
+    allocate(density(var%nz))
+
+    if (dat%z_file(dat%nzf) < var%z(var%nz)) then
+      ! If TOA file is smaller than TOA of model.
+      ! We must extrapolate density of the file to higher altitudes
+      z_file = [z_file, var%z(var%nz)]
+      slope = (log10(dat%den_file(dat%nzf)) - log10(dat%den_file(dat%nzf-1)))/(dat%z_file(dat%nzf) - dat%z_file(dat%nzf-1))
+      y_intercept = log10(dat%den_file(dat%nzf)) - slope*dat%z_file(dat%nzf)
+      den_tmp = slope*var%z(var%nz) + y_intercept
+      den_tmp = 10.0_dp**den_tmp
+      den_file = [den_file, den_tmp]
+      nzf1 = nzf1 + 1
+    elseif (dat%z_file(1) > var%z(1)) then
+      ! If BOA file is high altitude then BOA of model.
+      ! We must extrapolate density of the file to lower altitudes
+      z_file = [var%z(1), z_file]
+      slope = (log10(dat%den_file(2)) - log10(dat%den_file(1)))/(dat%z_file(2) - dat%z_file(1))
+      y_intercept = log10(dat%den_file(1)) - slope*dat%z_file(1)
+      den_tmp = slope*var%z(1) + y_intercept
+      den_tmp = 10.0_dp**den_tmp
+      den_file = [den_tmp, den_file]
+      nzf1 = nzf1 + 1
+    endif
+
+    ! Interpolate file density to model grid
+    call interp(var%nz, nzf1, var%z, z_file,&
+                log10(den_file), density, ierr)
+    if (ierr /= 0) then
+      err = 'Subroutine interp returned an error.'
+      return
+    endif
+    density = 10.0_dp**density
+
+    do i = 1,dat%nq
+      call interp(var%nz, dat%nzf, var%z, dat%z_file,&
+                  log10(abs(dat%mix_file(i,:))), var%usol_init(i,:), ierr)
+      if (ierr /= 0) then
+        err = 'Subroutine interp returned an error.'
+        return
+      endif
+    enddo
+    var%usol_init = 10.0_dp**var%usol_init
+
+    do i = 1,var%nz
+      var%usol_init(:,i) = var%usol_init(:,i)*density(i)
+    enddo 
+
+  end subroutine
+
   module subroutine compute_gibbs_energy(dat, var, err)
     use photochem_eqns, only: gibbs_energy_eval
 

src/input/photochem_input_read.f90

@@ -720,6 +720,9 @@ subroutine unpack_settings(infile, s, dat, var, err)
     endif
 
     ! default-gas-lower-boundary already applied to PhotoSettings
+
+    ! conserving initialization?
+    dat%conserving_init = s%conserving_init
 
     ! water
     dat%fix_water_in_trop = s%fix_water_in_trop

src/photochem_types.f90

@@ -64,6 +64,8 @@ module photochem_types ! make a giant IO object
     integer :: H_escape_type
     real(dp), allocatable :: H_escape_S1
     integer :: default_lowerboundcond
+    ! initialization
+    logical :: conserving_init
     ! climate
     logical :: evolve_climate
     ! water
@@ -334,6 +336,10 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
     type(ParticleXsections), allocatable :: part_xs_file(:) !! np in length
 
     ! initial conditions  
+    !> Only used in EVOATMOSPHERE.
+    !> If True, then during initialization, molecules in the atmosphere will
+    !> be conserved.
+    logical :: conserving_init = .false.
     integer :: nzf !! number of atmospheric layers in file
     real(dp), allocatable :: z_file(:) !! (nzf) cm
     real(dp), allocatable :: T_file(:) !! (nzf) K

src/photochem_types_create.f90

@@ -63,6 +63,11 @@ function unpack_PhotoSettings(root, filename, err) result(s)
     s%nz = dict%get_integer('number-of-layers',error = io_err)
     if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
 
+    if (s%bottom /= 0.0_dp) then
+      err = 'The bottom of the atmosphere must be 0.'
+      return
+    endif
+
     if (s%top < s%bottom) then
       err = 'The top of the atmosphere must be bigger than the bottom'
       return
@@ -193,6 +198,10 @@ function unpack_PhotoSettings(root, filename, err) result(s)
       return
     endif
 
+    ! Atmosphere initialization
+    s%conserving_init = dict%get_logical('conserving-initialization',default=.false.,error = io_err)
+    if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
+
     ! climate
     s%evolve_climate = dict%get_logical('evolve-climate',default=.false.,error = io_err)
     if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif