allow ability to initialize TOA from file

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION "3.14")
 cmake_policy(SET CMP0148 OLD)
-project(Photochem LANGUAGES Fortran C VERSION "0.5.2")
+project(Photochem LANGUAGES Fortran C VERSION "0.5.3")
 
 include(FortranCInterface)
 FortranCInterface_VERIFY()

src/input/photochem_input_read.f90

@@ -29,6 +29,10 @@ module subroutine read_all_files(mechanism_file, s, flux_file, atmosphere_txt, b
 
     call get_photomech(mechanism_file, dat, var, err)
     if (allocated(err)) return
+
+    ! initial atmosphere
+    call read_atmosphere_file(atmosphere_txt, dat, var, err)
+    if (allocated(err)) return
 
     call unpack_settings(s%filename, s, dat, var, err)
     if (allocated(err)) return
@@ -47,10 +51,6 @@ module subroutine read_all_files(mechanism_file, s, flux_file, atmosphere_txt, b
     call read_stellar_flux(flux_file, dat%nw, dat%wavl, var%photon_flux, err)
     if (allocated(err)) return
 
-    ! initial atmosphere
-    call read_atmosphere_file(atmosphere_txt, dat, var, err)
-    if (allocated(err)) return
-
   end subroutine
 
   subroutine get_photomech(infile, dat, var, err) 
@@ -639,7 +639,22 @@ subroutine unpack_settings(infile, s, dat, var, err)
     !!! atmosphere-grid !!!
     !!!!!!!!!!!!!!!!!!!!!!!
     var%bottom_atmos = s%bottom
-    var%top_atmos = s%top
+
+    read(s%top,*,iostat = io) var%top_atmos
+    if (io /= 0) then
+      ! it isn't a float
+      if (trim(s%top) == "atmospherefile") then
+        ! Use the TOA in the initialization file
+        var%top_atmos = dat%z_file(dat%nzf) + 0.5_dp*(dat%z_file(2) - dat%z_file(1))
+      else
+        err = '"top" can only be a number or "atmospherefile". See '//trim(infile)
+        return 
+      endif
+    endif
+    if (var%top_atmos < var%bottom_atmos) then
+      err = 'The top of the atmosphere must be bigger than the bottom'
+      return
+    endif
     var%nz = s%nz
 
     !!!!!!!!!!!!!!!!!!!!!!!
@@ -772,6 +787,10 @@ subroutine unpack_settings(infile, s, dat, var, err)
     if (dat%fix_water_in_trop .or. dat%gas_rainout) then
       ! we need a tropopause altitude
       var%trop_alt = s%trop_alt
+      if (var%trop_alt > var%top_atmos) then
+        err = 'IOError: tropopause-altitude must be between the top and bottom of the atmosphere'
+        return
+      endif
     endif
 
     if (dat%water_cond) then        
@@ -2532,7 +2551,6 @@ subroutine read_stellar_flux(star_file, nw, wavl, photon_flux, err)
 
   end subroutine
 
-
   subroutine read_atmosphere_file(atmosphere_txt, dat, var, err)
     character(len=*), intent(in) :: atmosphere_txt
     type(PhotochemData), intent(inout) :: dat

src/photochem_types.f90

@@ -43,7 +43,7 @@ module photochem_types ! make a giant IO object
 
     ! atmosphere-grid
     real(dp) :: bottom
-    real(dp) :: top
+    character(:), allocatable :: top
     integer :: nz
 
     ! photolysis-grid

src/photochem_types_create.f90

@@ -58,7 +58,7 @@ function unpack_PhotoSettings(root, filename, err) result(s)
     if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
     s%bottom = dict%get_real('bottom',error = io_err)
     if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
-    s%top = dict%get_real('top',error = io_err)
+    s%top = trim(dict%get_string('top',error = io_err))
     if (allocated(io_err)) then; err = trim(filename)//trim(io_err%message); return; endif
     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
@@ -68,11 +68,6 @@ function unpack_PhotoSettings(root, filename, err) result(s)
       return
     endif
 
-    if (s%top < s%bottom) then
-      err = 'The top of the atmosphere must be bigger than the bottom'
-      return
-    endif
-
     if (s%nz < 10) then
       err = "The number of vertical layers must be >= 10"
       return
@@ -256,10 +251,9 @@ function unpack_PhotoSettings(root, filename, err) result(s)
         err = "tropopause-altitude must be specified if fix-water-in-troposphere = true, or gas-rainout = true"
         return
       endif
-      if ((s%trop_alt < s%bottom) .or. &
-          (s%trop_alt > s%top)) then
-          err = 'IOError: tropopause-altitude must be between the top and bottom of the atmosphere'
-          return
+      if (s%trop_alt < s%bottom) then
+        err = 'IOError: tropopause-altitude must be between the top and bottom of the atmosphere'
+        return
       endif
 
     endif