From 7d9420c5e363a71e45965e3f540f345c71c33214 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Fri, 21 Jun 2024 09:05:56 -0700 Subject: [PATCH] faster RCE by only doing IR radiative transfer --- CMakeLists.txt | 2 +- src/adiabat/clima_adiabat.f90 | 3 +++ src/adiabat/clima_adiabat_rc.f90 | 4 ++-- src/adiabat/clima_adiabat_solve.f90 | 11 ++++++----- src/radtran/clima_radtran.f90 | 18 ++++++++++++++---- 5 files changed, 26 insertions(+), 12 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 658a358..99d8451 100644 --- a/CMakeLists.txt +++ b/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") diff --git a/src/adiabat/clima_adiabat.f90 b/src/adiabat/clima_adiabat.f90 index b7bbb26..e099d02 100644 --- a/src/adiabat/clima_adiabat.f90 +++ b/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 diff --git a/src/adiabat/clima_adiabat_rc.f90 b/src/adiabat/clima_adiabat_rc.f90 index f6fdce2..92ba542 100644 --- a/src/adiabat/clima_adiabat_rc.f90 +++ b/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 diff --git a/src/adiabat/clima_adiabat_solve.f90 b/src/adiabat/clima_adiabat_solve.f90 index e70ad25..b54d4d2 100644 --- a/src/adiabat/clima_adiabat_solve.f90 +++ b/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 diff --git a/src/radtran/clima_radtran.f90 b/src/radtran/clima_radtran.f90 index b790784..d409dea 100644 --- a/src/radtran/clima_radtran.f90 +++ b/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