faster RCE by only doing IR radiative transfer

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION "3.14")
 cmake_policy(SET CMP0148 OLD) # To suppress a warning emerging from scikit-build
 
-project(Clima LANGUAGES Fortran C VERSION "0.4.6")
+project(Clima LANGUAGES Fortran C VERSION "0.4.7")
 
 set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/modules")
 

src/adiabat/clima_adiabat.f90

@@ -107,6 +107,9 @@ module clima_adiabat
     !> A term that weights the importance of maintaining radiative
     !> equilibrium to convection.
     real(dp) :: radiation_norm_term = 1.0e-3_dp
+    !> If False, then the jacobian calculation in RCE does not recompute
+    !> solar radiative transfer.
+    logical :: compute_solar_in_jac = .false.
     logical :: verbose = .true. !! verbosity
 
     ! State of the atmosphere

src/adiabat/clima_adiabat_rc.f90

@@ -444,8 +444,8 @@ subroutine integrate(d, err)
       enddo
 
       if (k > 2*d%nz) then
-        err = 'Terrible out-of-bounds error in clima_adiabat_rc. Report this bug!'
-        return
+        print*,'Terrible out-of-bounds error in clima_adiabat_rc. Report this bug!'
+        stop 1
       endif
 
       if (k == 2*d%nz) exit

src/adiabat/clima_adiabat_solve.f90

@@ -281,17 +281,17 @@ subroutine AdiabatClimate_objective(self, P_i_surf, x, include_heat_capacity, re
 
     ! resets self%T_surf, self%T, self%densities, self%lapse_rate
     ! also does radiative transfer and computes res
-    call AdiabatClimate_objective_(self, P_i_surf, T_in, include_heat_capacity, res, err)
+    call AdiabatClimate_objective_(self, P_i_surf, T_in, include_heat_capacity, .true., res, err)
     if (allocated(err)) return
 
   end subroutine
 
-  subroutine AdiabatClimate_objective_(self, P_i_surf, T_in, include_heat_capacity, res, err)
+  subroutine AdiabatClimate_objective_(self, P_i_surf, T_in, include_heat_capacity, compute_solar, res, err)
     use clima_const, only: k_boltz
     class(AdiabatClimate), intent(inout) :: self
     real(dp), intent(in) :: P_i_surf(:)
     real(dp), intent(in) :: T_in(:)
-    logical, intent(in) :: include_heat_capacity
+    logical, intent(in) :: include_heat_capacity, compute_solar
     real(dp), intent(out) :: res(:)
     character(:), allocatable, intent(out) :: err
 
@@ -323,7 +323,8 @@ subroutine AdiabatClimate_objective_(self, P_i_surf, T_in, include_heat_capacity
 
     ! radiative transfer
     call self%copy_atm_to_radiative_grid()
-    call self%rad%radiate(self%T_surf, self%T_r, self%P_r/1.0e6_dp, self%densities_r, self%dz_r, err=err)
+    call self%rad%radiate(self%T_surf, self%T_r, self%P_r/1.0e6_dp, self%densities_r, self%dz_r, &
+      compute_solar=compute_solar, err=err)
     if (allocated(err)) return
     if (self%tidally_locked_dayside) then; block
       real(dp) :: tau_LW, k_term, f_term, rad_enhancement
@@ -393,7 +394,7 @@ subroutine AdiabatClimate_jacobian(self, P_i_surf, x, include_heat_capacity, jac
         T_in(self%ind_conv_lower(ind):self%ind_conv_upper(ind)) + deltaT
       endif
 
-      call AdiabatClimate_objective_(self, P_i_surf, T_perturb, include_heat_capacity, res_perturb, err)
+      call AdiabatClimate_objective_(self, P_i_surf, T_perturb, include_heat_capacity, self%compute_solar_in_jac, res_perturb, err)
       if (allocated(err)) return
 
       ! Compute jacobian

src/radtran/clima_radtran.f90

@@ -199,7 +199,7 @@ function create_Radtran_2(species_names, particle_names, s, star_f, &
 
   end function
 
-  subroutine Radtran_radiate(self, T_surface, T, P, densities, dz, pdensities, radii , err)
+  subroutine Radtran_radiate(self, T_surface, T, P, densities, dz, pdensities, radii, compute_solar, err)
     use clima_radtran_radiate, only: radiate
     use clima_const, only: pi
     class(Radtran), target, intent(inout) :: self
@@ -209,20 +209,27 @@ subroutine Radtran_radiate(self, T_surface, T, P, densities, dz, pdensities, rad
     real(dp), intent(in) :: densities(:,:) !! (nz,ng) number density of each 
                                            !! molecule in each layer (molcules/cm3)
     real(dp), optional, target, intent(in) :: pdensities(:,:), radii(:,:)
+    logical, optional, intent(in) :: compute_solar
     real(dp), intent(in) :: dz(:) !! (nz) thickness of each layer (cm)
     character(:), allocatable, intent(out) :: err
 
     integer :: ierr
+    logical :: compute_solar_
 
     type(ClimaRadtranWrk), pointer :: wrk_ir
     type(ClimaRadtranWrk), pointer :: wrk_sol
 
     wrk_ir => self%wrk_ir
-    wrk_sol => self%wrk_sol                                        
+    wrk_sol => self%wrk_sol                                    
 
     call check_inputs(self%nz, self%ng, self%np, T, P, densities, dz, pdensities, radii, err)
     if (allocated(err)) return
 
+    compute_solar_ = .true.
+    if (present(compute_solar)) then
+      compute_solar_ = compute_solar
+    endif
+
     ! IR radiative transfer                                     
     ierr = radiate(self%ir, &
                    self%surface_emissivity, &
@@ -236,6 +243,7 @@ subroutine Radtran_radiate(self, T_surface, T, P, densities, dz, pdensities, rad
       err = 'Input particle radii are outside the data range.'
       return
     endif
+    if (compute_solar_) then
     ! Solar radiative transfer
     ierr = radiate(self%sol, &
                    self%surface_emissivity, &
@@ -249,13 +257,14 @@ subroutine Radtran_radiate(self, T_surface, T, P, densities, dz, pdensities, rad
       err = 'Input particle radii are outside the data range.'
       return
     endif
+    endif
     ! Total flux at edges of layers (ergs/(cm2 s) which is the same as mW/m2).
     ! Index 1 is bottom. Index nz+1 is top edge of top layer.
     self%f_total = (wrk_sol%fdn_n - wrk_sol%fup_n) + (wrk_ir%fdn_n - wrk_ir%fup_n)
 
   end subroutine
 
-  subroutine Radtran_TOA_fluxes(self, T_surface, T, P, densities, dz, pdensities, radii, ISR, OLR, err)
+  subroutine Radtran_TOA_fluxes(self, T_surface, T, P, densities, dz, pdensities, radii, compute_solar, ISR, OLR, err)
     use clima_radtran_radiate, only: radiate
     class(Radtran), target, intent(inout) :: self
     real(dp), intent(in) :: T_surface
@@ -265,10 +274,11 @@ subroutine Radtran_TOA_fluxes(self, T_surface, T, P, densities, dz, pdensities,
                                            !! molecule in each layer (molcules/cm3)
     real(dp), intent(in) :: dz(:) !! (nz) thickness of each layer (cm)
     real(dp), optional, target, intent(in) :: pdensities(:,:), radii(:,:) !! (nz,np)
+    logical, optional, intent(in) :: compute_solar
     real(dp), intent(out) :: ISR, OLR
     character(:), allocatable, intent(out) :: err
 
-    call self%radiate(T_surface, T, P, densities, dz, pdensities, radii , err)
+    call self%radiate(T_surface, T, P, densities, dz, pdensities, radii, compute_solar, err)
     if (allocated(err)) return
 
     ISR = (self%wrk_sol%fdn_n(self%nz+1) - self%wrk_sol%fup_n(self%nz+1)) ! Incoming short wave