added ability to do tidally locked daysides

src/adiabat/clima_adiabat.f90

@@ -63,6 +63,16 @@ module clima_adiabat
     !> reservoir of gas dissolved in oceans in mol/cm^2 (ng, ng). There can be multiple oceans.
     !> The gases dissolved in ocean made of species 1 is given by `N_ocean(:,1)`.
     real(dp), allocatable :: N_ocean(:,:) 
+
+    ! Heat re-distribution terms for Equation (10) in Koll (2020), ApJ when
+    ! considering tidally locked planets. Default values are from the paper.
+    !> If true, then will attempt to compute the climate corresponding to the
+    !> observed dayside temperature of a tidally locked planet.
+    logical :: tidally_locked_dayside = .false.
+    real(dp) :: L !! = planet radius. Circulation’s horizontal scale (cm)
+    real(dp) :: chi = 0.2_dp !! Heat engine efficiency term (no units)
+    real(dp) :: n_LW = 2.0_dp !! = 1 or 2 (no units)
+    real(dp) :: Cd = 1.9e-3_dp !! Drag coefficient (no units)
 
   contains
     ! Constructs atmospheres
@@ -81,6 +91,8 @@ module clima_adiabat
     procedure :: set_ocean_solubility_fcn => AdiabatClimate_set_ocean_solubility_fcn
     procedure :: to_regular_grid => AdiabatClimate_to_regular_grid
     procedure :: out2atmosphere_txt => AdiabatClimate_out2atmosphere_txt
+    ! For tidally locked planets
+    procedure :: heat_redistribution_parameters => AdiabatClimate_heat_redistribution_parameters
   end type
 
   interface AdiabatClimate
@@ -164,6 +176,9 @@ function create_AdiabatClimate(species_f, settings_f, star_f, data_dir, err) res
     allocate(c%P(c%nz), c%T(c%nz), c%f_i(c%nz,c%sp%ng), c%z(c%nz), c%dz(c%nz))
     allocate(c%densities(c%nz,c%sp%ng))
     allocate(c%N_atmos(c%sp%ng),c%N_surface(c%sp%ng),c%N_ocean(c%sp%ng,c%sp%ng))
+
+    ! Heat redistribution parameter
+    c%L = c%planet_radius
 
   end function
 
@@ -472,7 +487,7 @@ subroutine fcn(n_, x_, fvec_, iflag_)
       real(dp), intent(in) :: x_(n_)
       real(dp), intent(out) :: fvec_(n_)
       integer, intent(inout) :: iflag_
-      real(dp) :: T, ISR, OLR
+      real(dp) :: T, ISR, OLR, rad_enhancement
       T = 10.0_dp**x_(1)
       if (self%solve_for_T_trop) then
         self%T_trop = 10.0_dp**x_(2)
@@ -482,12 +497,31 @@ subroutine fcn(n_, x_, fvec_, iflag_)
         iflag_ = -1
         return
       endif
-      fvec_(1) = ISR - OLR
+      rad_enhancement = 1.0_dp
+      if (self%tidally_locked_dayside) then; block
+        real(dp) :: tau_LW, k_term, f_term
+        call self%heat_redistribution_parameters(tau_LW, k_term, f_term, err)
+        if (allocated(err)) then
+          iflag_ = -1
+          return
+        endif
+        ! Increase the stellar flux, because we are computing the climate of
+        ! observed dayside.
+        rad_enhancement = 4.0_dp*f_term
+      endblock; endif
+      fvec_(1) = ISR*rad_enhancement - OLR
 
       if (self%solve_for_T_trop) then; block
-        real(dp) :: bond_albedo
+        use clima_eqns, only: skin_temperature
+        real(dp) :: bond_albedo, stellar_radiation
+        integer :: i
         bond_albedo = self%rad%wrk_sol%fup_n(self%nz+1)/self%rad%wrk_sol%fdn_n(self%nz+1)
-        fvec_(2) = self%rad%skin_temperature(bond_albedo) - self%T_trop
+        stellar_radiation = 0.0_dp
+        do i = 1,self%rad%sol%nw
+          stellar_radiation = stellar_radiation + self%rad%photons_sol(i)*(self%rad%sol%freq(i) - self%rad%sol%freq(i+1))
+        enddo
+        stellar_radiation = stellar_radiation/1.0e3_dp
+        fvec_(2) = skin_temperature(stellar_radiation*rad_enhancement, bond_albedo) - self%T_trop
       endblock; endif
     end subroutine
 
@@ -538,7 +572,7 @@ subroutine fcn(n_, x_, fvec_, iflag_)
       real(dp), intent(in) :: x_(n_)
       real(dp), intent(out) :: fvec_(n_)
       integer, intent(inout) :: iflag_
-      real(dp) :: T, ISR, OLR
+      real(dp) :: T, ISR, OLR, rad_enhancement
       T = 10.0_dp**x_(1)
       if (self%solve_for_T_trop) then
         self%T_trop = 10.0_dp**x_(2)
@@ -548,12 +582,31 @@ subroutine fcn(n_, x_, fvec_, iflag_)
         iflag_ = -1
         return
       endif
-      fvec_(1) = ISR - OLR
+      rad_enhancement = 1.0_dp
+      if (self%tidally_locked_dayside) then; block
+        real(dp) :: tau_LW, k_term, f_term
+        call self%heat_redistribution_parameters(tau_LW, k_term, f_term, err)
+        if (allocated(err)) then
+          iflag_ = -1
+          return
+        endif
+        ! Increase the stellar flux, because we are computing the climate of
+        ! observed dayside.
+        rad_enhancement = 4.0_dp*f_term
+      endblock; endif
+      fvec_(1) = ISR*rad_enhancement - OLR
 
       if (self%solve_for_T_trop) then; block
-        real(dp) :: bond_albedo
+        use clima_eqns, only: skin_temperature
+        real(dp) :: bond_albedo, stellar_radiation
+        integer :: i
         bond_albedo = self%rad%wrk_sol%fup_n(self%nz+1)/self%rad%wrk_sol%fdn_n(self%nz+1)
-        fvec_(2) = self%rad%skin_temperature(bond_albedo) - self%T_trop
+        stellar_radiation = 0.0_dp
+        do i = 1,self%rad%sol%nw
+          stellar_radiation = stellar_radiation + self%rad%photons_sol(i)*(self%rad%sol%freq(i) - self%rad%sol%freq(i+1))
+        enddo
+        stellar_radiation = stellar_radiation/1.0e3_dp
+        fvec_(2) = skin_temperature(stellar_radiation*rad_enhancement, bond_albedo) - self%T_trop
       endblock; endif
     end subroutine
 
@@ -607,7 +660,7 @@ subroutine fcn(n_, x_, fvec_, iflag_)
       real(dp), intent(in) :: x_(n_)
       real(dp), intent(out) :: fvec_(n_)
       integer, intent(inout) :: iflag_
-      real(dp) :: T, ISR, OLR
+      real(dp) :: T, ISR, OLR, rad_enhancement
       T = 10.0_dp**x_(1)
       if (self%solve_for_T_trop) then
         self%T_trop = 10.0_dp**x_(2)
@@ -617,12 +670,31 @@ subroutine fcn(n_, x_, fvec_, iflag_)
         iflag_ = -1
         return
       endif
-      fvec_(1) = ISR - OLR
+      rad_enhancement = 1.0_dp
+      if (self%tidally_locked_dayside) then; block
+        real(dp) :: tau_LW, k_term, f_term
+        call self%heat_redistribution_parameters(tau_LW, k_term, f_term, err)
+        if (allocated(err)) then
+          iflag_ = -1
+          return
+        endif
+        ! Increase the stellar flux, because we are computing the climate of
+        ! observed dayside.
+        rad_enhancement = 4.0_dp*f_term
+      endblock; endif
+      fvec_(1) = ISR*rad_enhancement - OLR
 
       if (self%solve_for_T_trop) then; block
-        real(dp) :: bond_albedo
+        use clima_eqns, only: skin_temperature
+        real(dp) :: bond_albedo, stellar_radiation
+        integer :: i
         bond_albedo = self%rad%wrk_sol%fup_n(self%nz+1)/self%rad%wrk_sol%fdn_n(self%nz+1)
-        fvec_(2) = self%rad%skin_temperature(bond_albedo) - self%T_trop
+        stellar_radiation = 0.0_dp
+        do i = 1,self%rad%sol%nw
+          stellar_radiation = stellar_radiation + self%rad%photons_sol(i)*(self%rad%sol%freq(i) - self%rad%sol%freq(i+1))
+        enddo
+        stellar_radiation = stellar_radiation/1.0e3_dp
+        fvec_(2) = skin_temperature(stellar_radiation*rad_enhancement, bond_albedo) - self%T_trop
       endblock; endif
     end subroutine
   end function
@@ -784,5 +856,83 @@ subroutine AdiabatClimate_out2atmosphere_txt(self, filename, eddy, overwrite, cl
     close(2)
 
   end subroutine
+
+  !> For considering a tidally locked planet. This function computes key parameters for
+  !> Equation (10) in Koll (2022), ApJ. The function must be called after calling a function
+  !> `self%TOA_fluxes`, because it uses atmosphere properties, and radiative properties.
+  subroutine AdiabatClimate_heat_redistribution_parameters(self, tau_LW, k_term, f_term, err)
+    use clima_const, only: c_light
+    use clima_eqns, only: k_term_heat_redistribution, f_heat_redistribution, &
+                          gravity, heat_capacity_eval, planck_fcn
+    class(AdiabatClimate), intent(inout) :: self
+    real(dp), intent(out) :: tau_LW
+    real(dp), intent(out) :: k_term
+    real(dp), intent(out) :: f_term
+    character(:), allocatable, intent(out) :: err
+
+    integer :: i
+
+    real(dp) :: bond_albedo, Teq
+    real(dp) :: grav
+    real(dp) :: mubar
+    logical :: found
+    real(dp) ::  cp_tmp, cp
+
+    ! equilibrium temperature
+    bond_albedo = self%rad%wrk_sol%fup_n(self%nz+1)/self%rad%wrk_sol%fdn_n(self%nz+1)
+    Teq = self%rad%equilibrium_temperature(bond_albedo)
+
+    ! gravity
+    grav = gravity(self%planet_radius, self%planet_mass, 0.0_dp)
+
+    ! mean molecular weight at surface
+    mubar = 0.0_dp
+    do i = 1,self%sp%ng
+      mubar = mubar + self%f_i(1,i)*self%sp%g(i)%mass
+    enddo
+
+    ! heat capacity at surface
+    cp = tiny(0.0_dp)
+    do i = 1,self%sp%ng
+      call heat_capacity_eval(self%sp%g(i)%thermo, self%T_surf, found, cp_tmp) ! J/(mol*K)
+      if (.not. found) then
+        err = "Failed to compute heat capacity"
+        return
+      endif
+      ! J/(mol*K)
+      cp = cp + cp_tmp*self%f_i(1,i) ! J/(mol*K)
+    enddo
+    ! J/(mol*K) * (mol/kg) = J/(kg*K)
+    cp = cp*(1.0_dp/(mubar*1.0e-3_dp))
+    ! J/(kg*K) * (erg/J) * (kg/g) = erg/(g*K)
+    cp = cp*1.0e4_dp
+
+    ! tau_LW. Computed with Equation (13) in Koll (2020), ApJ
+    block
+      real(dp) :: numerator, denominator, dlambda, tau_lambda, bplank, avg_freq, avg_lam
+      numerator = 0.0_dp
+      denominator = 0.0_dp
+      do i = 1,self%rad%ir%nw
+        dlambda = self%rad%ir%wavl(i+1) - self%rad%ir%wavl(i) ! Width of a wavelength bin (nm)
+        tau_lambda = sum(self%rad%wrk_ir%rx%tau_band(:,i)) ! IR optical depth
+
+        avg_freq = 0.5_dp*(self%rad%ir%freq(i) + self%rad%ir%freq(i+1))
+        avg_lam = (c_light*1.0e9_dp/avg_freq) ! (m/s) * (nm/m) * (s) = nm
+
+        bplank = planck_fcn(avg_freq, self%T_surf) ! (mW sr^−1 m^−2 Hz^-1)
+        bplank = bplank * (avg_freq/avg_lam) ! (mW sr^−1 m^−2 Hz^-1) * (Hz/nm) = (mW sr^−1 m^−2 nm^-1)
+
+        ! integrate the numerator
+        numerator = numerator + exp(-tau_lambda)*bplank*dlambda
+        ! integrate the denominator
+        denominator = denominator + bplank*dlambda
+      enddo
+      tau_LW = - log(numerator/denominator)
+    endblock
+
+    k_term = k_term_heat_redistribution(self%L, grav, self%chi, mubar, cp, self%n_LW, self%Cd)
+    f_term = f_heat_redistribution(tau_LW, self%P_surf, Teq, k_term)
+
+  end subroutine
 
 end module

src/clima_eqns.f90

@@ -250,4 +250,41 @@ pure function skin_temperature(stellar_radiation, bond_albedo) result(T_skin)
     T_skin = equilibrium_temperature(stellar_radiation, bond_albedo)*(0.5_dp)**(0.25_dp)
   end function
 
+  !> k term in Equation (10) of Koll (2022), ApJ
+  pure function k_term_heat_redistribution(L, grav, chi, mubar, cp, n_LW, Cd) result(k)
+    use clima_const, only: sigma_si, Rgas
+    real(dp), intent(in) :: L !! Circulation’s horizontal scale (cm)
+    real(dp), intent(in) :: grav !! Surface gravity (cm/s^2)
+    real(dp), intent(in) :: chi !! Heat engine efficiency term (no units)
+    real(dp), intent(in) :: mubar !! Mean molar weight (g/mol)
+    real(dp), intent(in) :: cp !! Heat capacity (erg/(g*K))
+    real(dp), intent(in) :: n_LW !! = 1 or 2 (no units)
+    real(dp), intent(in) :: Cd !! Drag coefficient (no units)
+    real(dp) :: k !! k term
+
+    real(dp) :: R_bar, Beta
+    real(dp), parameter :: sigma_cgs = sigma_si*1.0e3_dp ! Convert to cgs
+
+    R_bar = Rgas/mubar
+    Beta = R_bar/(cp*n_LW)
+    k = (L*grav)/(chi*Beta*cp) * ((Cd*sigma_cgs**2.0_dp)/R_bar)**(1.0_dp/3.0_dp) &
+        * (1.0e6_dp)**(-2.0_dp/3.0_dp) * (600.0_dp)**(4.0_dp/3.0_dp)
+
+  end function
+
+  !> The heat redistribution parameter `f` from Equation (10) in Koll (2022), ApJ.
+  pure function f_heat_redistribution(tau_LW, Ps, Teq, k) result(f)
+    real(dp), intent(in) :: tau_LW !! Long wavelength optical depth
+    real(dp), intent(in) :: Ps !! Surface pressure (dynes/cm^2)
+    real(dp), intent(in) :: Teq !! Planetary equilibrium temperature (K)
+    real(dp), intent(in) :: k !! k term
+
+    real(dp) :: f
+
+    f = (2.0_dp/3.0_dp) - (5.0_dp/12.0_dp) &
+        * (tau_LW**(1.0_dp/3.0_dp) * (Ps/1.0e6_dp)**(2.0_dp/3.0_dp)*(Teq/600.0_dp)**(-4.0_dp/3.0_dp)) &
+        / (k + tau_LW**(1.0_dp/3.0_dp) * (Ps/1.0e6_dp)**(2.0_dp/3.0_dp)*(Teq/600.0_dp)**(-4.0_dp/3.0_dp))
+
+  end function
+
 end module

src/radtran/clima_radtran.f90

@@ -89,6 +89,7 @@ module clima_radtran
     procedure :: radiate => Radtran_radiate
     procedure :: TOA_fluxes => Radtran_TOA_fluxes
     procedure :: skin_temperature => Radtran_skin_temperature
+    procedure :: equilibrium_temperature => Radtran_equilibrium_temperature
   end type
 
   interface Radtran
@@ -431,6 +432,23 @@ function Radtran_skin_temperature(self, bond_albedo) result(T_skin)
     T_skin = skin_temperature(stellar_radiation, bond_albedo)
   end function
 
+  function Radtran_equilibrium_temperature(self, bond_albedo) result(T_eq)
+    use clima_eqns, only: equilibrium_temperature
+    class(Radtran), target, intent(inout) :: self
+    real(dp), intent(in) :: bond_albedo
+    real(dp) :: T_eq
+
+    real(dp) :: stellar_radiation
+    integer :: i
+
+    stellar_radiation = 0.0_dp
+    do i = 1,self%sol%nw
+      stellar_radiation = stellar_radiation + self%photons_sol(i)*(self%sol%freq(i) - self%sol%freq(i+1))
+    enddo
+    stellar_radiation = stellar_radiation/1.0e3_dp
+    T_eq = equilibrium_temperature(stellar_radiation, bond_albedo)
+  end function
+
   !!!!!!!!!!!!!!!!!
   !!! Utilities !!!
   !!!!!!!!!!!!!!!!!

tests/test_adiabat.f90

@@ -21,6 +21,7 @@ program test
   endif
 
   c%solve_for_T_trop = .true.
+  c%tidally_locked_dayside = .true.
 
   T = c%surface_temperature( &
       [270.0e6_dp, 400e-6_dp*1.0e6_dp, 1.0e6_dp], &
@@ -31,6 +32,7 @@ program test
   endif
 
   print*,T, c%T_trop
+  c%tidally_locked_dayside = .false.
 
   T = c%surface_temperature_column( &
       [15.0e3_dp, 400e-6_dp*23.0_dp, 1.0*36.0_dp], &