From 07833d6ed149d118ed912e139db04d7af791b4cb Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Sat, 17 Feb 2024 16:01:18 -0800 Subject: [PATCH] make_profile_rc --- src/CMakeLists.txt | 1 + src/adiabat/clima_adiabat.f90 | 62 ++- src/adiabat/clima_adiabat_general.f90 | 1 + src/adiabat/clima_adiabat_rc.f90 | 666 ++++++++++++++++++++++++++ 4 files changed, 729 insertions(+), 1 deletion(-) create mode 100644 src/adiabat/clima_adiabat_rc.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5bb6f7c..ad86f09 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,6 +24,7 @@ set(RADIATE_SOURCES set(ADIABAT_SOURCES adiabat/clima_adiabat_general.f90 + adiabat/clima_adiabat_rc.f90 adiabat/clima_adiabat.f90 ) diff --git a/src/adiabat/clima_adiabat.f90 b/src/adiabat/clima_adiabat.f90 index 9f0100c..34d3d4d 100644 --- a/src/adiabat/clima_adiabat.f90 +++ b/src/adiabat/clima_adiabat.f90 @@ -14,7 +14,7 @@ module clima_adiabat ! settings and free parameters integer :: nz - real(dp) :: P_top = 1.0e-2_dp !! (dynes/cm2) + real(dp) :: P_top = 1.0_dp !! (dynes/cm2) real(dp) :: T_trop = 180.0_dp !! (T) real(dp), allocatable :: RH(:) !! relative humidity (ng) !> If .true., then any function that calls `make_column` will @@ -105,6 +105,8 @@ module clima_adiabat procedure :: surface_temperature => AdiabatClimate_surface_temperature procedure :: surface_temperature_column => AdiabatClimate_surface_temperature_column procedure :: surface_temperature_bg_gas => AdiabatClimate_surface_temperature_bg_gas + ! Test + procedure :: make_profile_rc => AdiabatClimate_make_profile_rc ! Utilities procedure :: set_ocean_solubility_fcn => AdiabatClimate_set_ocean_solubility_fcn procedure :: to_regular_grid => AdiabatClimate_to_regular_grid @@ -765,6 +767,64 @@ subroutine fcn(n_, x_, fvec_, iflag_) end subroutine end function + subroutine AdiabatClimate_make_profile_rc(self, P_i_surf, T_surf, T, err) + use clima_const, only: k_boltz, N_avo + use clima_adiabat_rc, only: make_profile_rc + class(AdiabatClimate), intent(inout) :: self + real(dp), intent(in) :: P_i_surf(:) !! dynes/cm^2 + real(dp), intent(in) :: T_surf, T(:) + character(:), allocatable, intent(out) :: err + + real(dp), allocatable :: P_e(:), z_e(:), f_i_e(:,:), lapse_rate_e(:) + real(dp), allocatable :: density(:) + integer :: i, j + + allocate(P_e(2*self%nz+1),z_e(2*self%nz+1),f_i_e(2*self%nz+1,self%sp%ng),lapse_rate_e(2*self%nz+1)) + allocate(density(self%nz)) + + if (size(P_i_surf) /= self%sp%ng) then + err = "P_i_surf has the wrong dimension" + return + endif + if (size(T) /= self%nz) then + err = "T has the wrong dimension" + endif + + self%T_surf = T_surf + self%T = T + call make_profile_rc(self%T_surf, self%T, P_i_surf, & + self%sp, self%nz, self%planet_mass, & + self%planet_radius, self%P_top, self%RH, & + self%rtol, self%atol, & + self%ocean_fcns, self%ocean_args_p, & + P_e, z_e, f_i_e, lapse_rate_e, & + self%N_surface, self%N_ocean, & + err) + if (allocated(err)) return + + self%T_surf = self%T_surf + self%P_surf = P_e(1) + do i = 1,self%nz + self%P(i) = P_e(2*i) + self%z(i) = z_e(2*i) + self%dz(i) = z_e(2*i+1) - z_e(2*i-1) + do j =1,self%sp%ng + self%f_i(i,j) = f_i_e(2*i,j) + enddo + enddo + + density = self%P/(k_boltz*self%T) + do j =1,self%sp%ng + self%densities(:,j) = self%f_i(:,j)*density(:) + enddo + + do i = 1,self%sp%ng + ! mol/cm^2 in atmosphere + self%N_atmos(i) = sum(density*self%f_i(:,i)*self%dz)/N_avo + enddo + + end subroutine + !> Sets a function for describing how gases dissolve in a liquid ocean. subroutine AdiabatClimate_set_ocean_solubility_fcn(self, species, fcn, err) class(AdiabatClimate), intent(inout) :: self diff --git a/src/adiabat/clima_adiabat_general.f90 b/src/adiabat/clima_adiabat_general.f90 index 2dcb3fe..c7de705 100644 --- a/src/adiabat/clima_adiabat_general.f90 +++ b/src/adiabat/clima_adiabat_general.f90 @@ -9,6 +9,7 @@ module clima_adiabat_general public :: make_profile, make_column public :: OceanFunction + public :: DrySpeciesType, CondensingSpeciesType type :: OceanFunction procedure(ocean_solubility_fcn), nopass, pointer :: fcn => null() diff --git a/src/adiabat/clima_adiabat_rc.f90 b/src/adiabat/clima_adiabat_rc.f90 new file mode 100644 index 0000000..4b1377b --- /dev/null +++ b/src/adiabat/clima_adiabat_rc.f90 @@ -0,0 +1,666 @@ +module clima_adiabat_rc + use clima_const, only: dp + use clima_types, only: Species + use dop853_module, only: dop853_class + use futils, only: brent_class + use clima_eqns, only: ocean_solubility_fcn + use clima_adiabat_general, only: OceanFunction, DrySpeciesType, CondensingSpeciesType + use linear_interpolation_module, only: linear_interp_1d + implicit none + + type :: AdiabatRCProfileData + ! Input data + real(dp), pointer :: T_surf + type(Species), pointer :: sp + integer, pointer :: nz + real(dp), pointer :: planet_mass + real(dp), pointer :: planet_radius + real(dp), pointer :: P_top + real(dp), pointer :: RH(:) + real(dp), pointer :: rtol + real(dp), pointer :: atol + ! Ouput + real(dp), pointer :: P(:) + real(dp), pointer :: z(:) + real(dp), pointer :: f_i(:,:) + real(dp), pointer :: lapse_rate(:) + + ! For interpolating the input T profile. + type(linear_interp_1d) :: T + real(dp), allocatable :: P_interp(:) + real(dp), allocatable :: T_interp(:) + + real(dp) :: P_surf !! surface pressure from inputs + + !> indicates whether species is dry or condensing (length ng) + integer, allocatable :: sp_type(:) + integer :: stopping_reason !! Reason why we stop integration + !> Index of root + integer :: ind_root + + !> Work space for root finding. All length (ng) + real(dp), allocatable :: gout(:) + real(dp), allocatable :: gout_old(:) + real(dp), allocatable :: gout_tmp(:) + logical, allocatable :: root_found(:) + real(dp), allocatable :: P_roots(:) + !> When we find a root, and exit integration, these + !> guys will give use the root + real(dp) :: P_root = huge(1.0_dp) + real(dp) :: u_root(1) + real(dp) :: epsj = 1.0e-8 + + !> Keeps track of how the dry atmosphere is partitioned + !> between different species + real(dp), allocatable :: f_i_dry(:) + + !> index for saving z, mixing ratios, and lapse rate + integer :: j + + ! work variables. All dimension (ng) + real(dp), allocatable :: P_i_cur(:) + real(dp), allocatable :: f_i_cur(:) + real(dp), allocatable :: cp_i_cur(:) + real(dp), allocatable :: L_i_cur(:) + + !> This helps us propogate error messages outside of the integration + character(:), allocatable :: err + + end type + + ! stopping_reason + enum, bind(c) + enumerator :: ReachedPtop, ReachedRoot + end enum + + type, extends(dop853_class) :: dop853_rc + type(AdiabatRCProfileData), pointer :: d => NULL() + end type + +contains + + ! inputs: P_i_surf, T_surf, T in each grid cell, P_top, planet_mass, planet_radius, nz, sp, ocean_fcns + ! indexes where convection is occuring + + ! outputs: Pressure, mixing ratios, z, lapse rate everywhere + + subroutine make_profile_rc(T_surf, T, P_i_surf, & + sp, nz, planet_mass, & + planet_radius, P_top, RH, & + rtol, atol, & + ocean_fcns, args_p, & + P, z, f_i, lapse_rate, & + N_surface, N_ocean, & + err) + use futils, only: linspace + use clima_eqns, only: gravity, ocean_solubility_fcn + use iso_c_binding, only: c_ptr + + real(dp), target, intent(in) :: T_surf !! K + real(dp), target, intent(in) :: T(:) !! (nz) + real(dp), intent(in) :: P_i_surf(:) !! (ng) dynes/cm2 + + type(Species), target, intent(inout) :: sp + integer, target, intent(in) :: nz + real(dp), target, intent(in) :: planet_mass, planet_radius + real(dp), target, intent(in) :: P_top, RH(:) + real(dp), target, intent(in) :: rtol, atol + type(OceanFunction), intent(in) :: ocean_fcns(:) + type(c_ptr), value, intent(in) :: args_p + + real(dp), target, intent(out) :: P(:), z(:) + real(dp), target, intent(out) :: f_i(:,:) + real(dp), target, intent(out) :: lapse_rate(:) + real(dp), target, intent(out) :: N_surface(:), N_ocean(:,:) + character(:), allocatable, intent(out) :: err + + type(AdiabatRCProfileData) :: d + integer :: i, j, ierr + real(dp) :: P_sat, grav, P_surface_inventory + + ! check inputs + if (size(T) /= nz) then + err = 'make_profile: Input "T" has the wrong shape' + return + endif + if (any(P_i_surf < 0.0_dp)) then + err = 'make_profile: Surface pressures (input "P_i_surf") must be positive' + return + endif + if (size(P_i_surf) /= sp%ng) then + err = 'make_profile: Input "P_i_surf" has the wrong shape' + return + endif + if (size(RH) /= sp%ng) then + err = 'make_profile: Input "RH" has the wrong dimension.' + return + endif + if (size(P) /= 2*nz+1) then + err = 'make_profile: Input "P" has the wrong shape' + return + endif + if (size(z) /= 2*nz+1) then + err = 'make_profile: Input "z" has the wrong shape' + return + endif + if (size(f_i, 1) /= 2*nz+1 .or. size(f_i, 2) /= sp%ng) then + err = 'make_profile: Input "f_i" has the wrong shape' + return + endif + if (size(lapse_rate) /= 2*nz+1) then + err = 'make_profile: Input "lapse_rate" has the wrong shape' + return + endif + if (size(N_surface) /= sp%ng) then + err = 'make_profile: Input "N_surface" has the wrong dimension.' + return + endif + if (size(N_ocean,1) /= sp%ng .or. size(N_ocean,2) /= sp%ng) then + err = 'make_profile: Input "N_ocean" has the wrong dimension.' + return + endif + + ! associate + ! inputs + d%T_surf => T_surf + d%sp => sp + d%nz => nz + d%planet_mass => planet_mass + d%planet_radius => planet_radius + d%P_top => P_top + d%RH => RH + d%rtol => rtol + d%atol => atol + ! outputs + d%P => P + d%z => z + d%f_i => f_i + d%lapse_rate => lapse_rate + + ! allocate work memory + allocate(d%sp_type(sp%ng)) + allocate(d%gout(sp%ng)) + allocate(d%gout_old(sp%ng)) + allocate(d%gout_tmp(sp%ng)) + allocate(d%P_roots(sp%ng)) + allocate(d%root_found(sp%ng)) + allocate(d%f_i_dry(sp%ng)) + allocate(d%P_i_cur(sp%ng)) + allocate(d%f_i_cur(sp%ng)) + allocate(d%cp_i_cur(sp%ng)) + allocate(d%L_i_cur(sp%ng)) + + ! gravity at the surface of planet + grav = gravity(planet_radius, planet_mass, 0.0_dp) + + do i = 1,d%sp%ng + ! compute saturation vapor pressures + P_sat = huge(1.0_dp) + if (allocated(d%sp%g(i)%sat)) then + P_sat = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T_surf) + endif + + ! determine if species are condensing, or not + if (P_i_surf(i) > P_sat) then + d%P_i_cur(i) = P_sat + ! the pressure of the ocean is everything not in the atmosphere + P_surface_inventory = P_i_surf(i) - P_sat + ! The surface mol/cm^2 can then be computed + ! from the "surface pressure" + N_surface(i) = P_surface_inventory/(sp%g(i)%mass*grav) + d%sp_type(i) = CondensingSpeciesType + else + d%P_i_cur(i) = P_i_surf(i) + N_surface(i) = 0.0_dp ! no surface reservoir if not condensing at the surface + d%sp_type(i) = DrySpeciesType + endif + enddo + + ! compute gases dissolved in oceans. This allows for multiple ocean. + do j = 1,sp%ng + if (associated(ocean_fcns(j)%fcn)) then; block + real(dp), allocatable :: m_i_cur(:) + + ! Compute mol/kg of each gas dissolved in the ocean + allocate(m_i_cur(sp%ng)) + call ocean_fcns(j)%fcn(T_surf, sp%ng, d%P_i_cur/1.0e6_dp, m_i_cur, args_p) + + ! Compute mol/cm^2 of each gas dissolved in ocean + N_ocean(j,j) = 0.0_dp ! ocean can not dissolve into itself + do i = 1,sp%ng + if (i /= j) then + N_ocean(i,j) = m_i_cur(i)*N_surface(j)*(sp%g(j)%mass/1.0e3_dp) + endif + enddo + endblock; else + ! Nothing dissolved in ocean + N_ocean(:,j) = 0.0_dp + endif + enddo + + d%P_surf = sum(d%P_i_cur) ! total pressure + if (P_top > d%P_surf) then + err = 'make_profile: "P_top" is bigger than the surface pressure' + return + endif + d%f_i_cur(:) = d%P_i_cur(:)/d%P_surf ! mixing ratios + call update_f_i_dry(d, d%P_surf, d%f_i_cur) + + ! Make P profile + call linspace(log10(d%P_surf),log10(P_top),P) + P(:) = 10.0_dp**P(:) + P(1) = d%P_surf + P(2*nz+1) = P_top + + ! Initialize the T interpolator + allocate(d%T_interp(nz+1),d%P_interp(nz+1)) + d%T_interp(1) = T_surf + d%P_interp(1) = log10(d%P_surf) + do i = 1,nz + d%T_interp(i+1) = T(i) + d%P_interp(i+1) = log10(P(2*i)) + enddo + d%T_interp = d%T_interp(size(d%T_interp):1:-1) + d%P_interp = d%P_interp(size(d%P_interp):1:-1) + call d%T%initialize(d%P_interp, d%T_interp, ierr) + if (ierr /= 0) then + err = 'Failed to initialize interpolator in "make_profile_rc"' + return + endif + + call integrate(d, err) + if (allocated(err)) return + + end subroutine + + subroutine integrate(d, err) + type(AdiabatRCProfileData), target, intent(inout) :: d + character(:), allocatable, intent(out) :: err + + type(dop853_rc) :: dop + logical :: status_ok + integer :: idid + real(dp) :: Pn, u(1), T_root + character(6) :: tmp_char + + call dop%initialize(fcn=right_hand_side_dop, solout=solout_dop, n=1, & + iprint=0, icomp=[1], status_ok=status_ok) + dop%d => d + d%j = 2 + ! Surface values + d%lapse_rate(1) = general_adiabat_lapse_rate(d, d%P_surf) + d%f_i(1,:) = d%f_i_cur(:) + d%z(1) = 0.0_dp + d%stopping_reason = ReachedPtop + if (.not. status_ok) then + err = 'dop853 initialization failed' + return + endif + + ! intial conditions + u = [d%z(1)] + Pn = d%P_surf + do + call dop%integrate(Pn, u, d%P_top, [d%rtol], [d%atol], iout=2, idid=idid) + if (allocated(d%err)) then + err = d%err + return + endif + if (idid < 0) then + write(tmp_char,'(i6)') idid + err = 'dop853 integration failed: '//trim(tmp_char) + return + endif + + if (d%stopping_reason == ReachedPtop) then + exit + elseif (d%stopping_reason == ReachedRoot) then; block + real(dp) :: f_dry + ! A root was hit. We restart integration + Pn = d%P_root ! root pressure + u = d%u_root ! root altitude + call d%T%evaluate(log10(Pn), T_root) + call mixing_ratios(d, Pn, T_root, d%f_i_cur, f_dry) ! get mixing ratios at the root + if (d%sp_type(d%ind_root) == DrySpeciesType) then + ! A gas has reached saturation. + d%sp_type(d%ind_root) = CondensingSpeciesType + elseif (d%sp_type(d%ind_root) == CondensingSpeciesType) then + ! A gas has reached a cold trap + d%sp_type(:) = DrySpeciesType + endif + call update_f_i_dry(d, Pn, d%f_i_cur) + d%stopping_reason = ReachedPtop + + endblock; endif + + enddo + + end subroutine + + subroutine solout_dop(self, nr, xold, x, y, irtrn, xout) + class(dop853_class),intent(inout) :: self + integer,intent(in) :: nr + real(dp),intent(in) :: xold + real(dp),intent(in) :: x + real(dp),dimension(:),intent(in) :: y + integer,intent(inout) :: irtrn + real(dp),intent(out) :: xout + + type(AdiabatRCProfileData), pointer :: d + real(dp) :: z_cur, P_cur, P_old + real(dp) :: T_i, PP, f_dry + integer :: i + + P_old = xold + P_cur = x + z_cur = y(1) + + select type (self) + class is (dop853_rc) + d => self%d + end select + + if (allocated(d%err)) then + ! if there is an error (from RHS function) + ! we return immediately + irtrn = -10 + return + endif + + call root_fcn(d, P_cur, d%gout) + if (nr == 1) then + d%gout_old(:) = d%gout(:) + endif + do i = 1,size(d%gout_old) + if ((d%gout_old(i) < 0 .and. d%gout(i) > 0) .or. & + (d%gout_old(i) > 0 .and. d%gout(i) < 0)) then + d%root_found(i) = .true. + call find_root(self, d, P_old, P_cur, i, d%P_roots(i)) + if (allocated(d%err)) then + irtrn = -10 + return + endif + else + d%root_found(i) = .false. + d%P_roots(i) = 0.0_dp + endif + enddo + + if (any(d%root_found)) then + ! lets use the highest pressure root. + ! this makes sure we don't skip an interesting root + irtrn = -1 + d%ind_root = maxloc(d%P_roots,1) ! highest pressure root + d%P_root = d%P_roots(d%ind_root) + d%u_root(1) = self%contd8(1, d%P_root) + d%stopping_reason = ReachedRoot + endif + + d%gout_old(:) = d%gout(:) ! save for next step + + ! save the results + if (d%j <= size(d%P)) then + if (irtrn == -1) then + PP = d%P_root + else + PP = P_cur + endif + + if (d%P(d%j) <= P_old .and. d%P(d%j) >= PP) then + do while (d%P(d%j) >= PP) + d%z(d%j) = self%contd8(1, d%P(d%j)) + call d%T%evaluate(log10(d%P(d%j)), T_i) + call mixing_ratios(d, d%P(d%j), T_i, d%f_i(d%j,:), f_dry) + d%lapse_rate(d%j) = general_adiabat_lapse_rate(d, d%P(d%j)) + d%j = d%j + 1 + if (d%j > size(d%P)) exit + enddo + endif + + endif + + end subroutine + + subroutine find_root(dop, d, P_old, P_cur, ind, P_root) + use futils, only: brent_class + type(dop853_class), intent(inout) :: dop + type(AdiabatRCProfileData), intent(inout) :: d + real(dp), intent(in) :: P_old + real(dp), intent(in) :: P_cur + integer, intent(in) :: ind + real(dp), intent(out) :: P_root + + real(dp), parameter :: tol = 1.0e-8_dp + real(dp) :: fzero + integer :: info + + type(brent_class) :: b + + call b%set_function(fcn) + call b%find_zero(P_old, P_cur, tol, P_root, fzero, info) + if (info /= 0) then + d%err = 'brent failed in "find_root"' + return + endif + + contains + function fcn(me_, x_) result(f_) + class(brent_class), intent(inout) :: me_ + real(dp), intent(in) :: x_ + real(dp) :: f_ + call root_fcn(d, x_, d%gout_tmp) + f_ = d%gout_tmp(ind) + end function + end subroutine + + subroutine root_fcn(d, P, gout) + type(AdiabatRCProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp), intent(out) :: gout(:) + + real(dp) :: T, f_dry, P_sat, dTdP + integer :: i + + call d%T%evaluate(log10(P), T) + call mixing_ratios(d, P, T, d%f_i_cur, f_dry) + d%P_i_cur = d%f_i_cur*P + + if (any(d%sp_type == CondensingSpeciesType)) then + dTdP = dT_dP(d, P) + endif + + do i = 1,d%sp%ng + ! compute saturation vapor pressures + P_sat = huge(1.0_dp) + if (allocated(d%sp%g(i)%sat)) then + P_sat = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T) + endif + + if (d%sp_type(i) == CondensingSpeciesType) then + ! Search for cold trap (place where Temperature decreases with altitude + gout(i) = dTdP - 1.0e-8_dp + elseif (d%sp_type(i) == DrySpeciesType) then + ! Dry species can become condensing species if + ! they reach saturation. + gout(i) = d%P_i_cur(i)/P_sat - (1.0_dp + 1.0e-8_dp) + endif + + enddo + + end subroutine + + subroutine right_hand_side_dop(self, P, u, du) + use clima_eqns, only: gravity, heat_capacity_eval + class(dop853_class), intent(inout) :: self + real(dp), intent(in) :: P + real(dp), intent(in), dimension(:) :: u + real(dp), intent(out), dimension(:) :: du + + select type (self) + class is (dop853_rc) + call right_hand_side(self%d, P, u, du) + end select + + end subroutine + + subroutine update_f_i_dry(d, P, f_i_layer) + type(AdiabatRCProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp), intent(in) :: f_i_layer(:) + + integer :: i + real(dp) :: P_dry + + d%P_i_cur(:) = f_i_layer(:)*P + P_dry = 0.0_dp + do i = 1,d%sp%ng + if (d%sp_type(i) == DrySpeciesType) then + P_dry = P_dry + d%P_i_cur(i) + endif + enddo + d%f_i_dry(:) = d%P_i_cur(:)/P_dry + + end subroutine + + subroutine mixing_ratios(d, P, T, f_i_layer, f_dry) + type(AdiabatRCProfileData), intent(inout) :: d + real(dp), intent(in) :: P, T + real(dp), intent(out) :: f_i_layer(:), f_dry + + real(dp) :: f_moist + integer :: i + + ! moist mixing ratios + f_moist = 0.0_dp + do i = 1,d%sp%ng + if (d%sp_type(i) == CondensingSpeciesType) then + f_i_layer(i) = min(d%RH(i)*d%sp%g(i)%sat%sat_pressure(T)/P,1.0_dp) + f_moist = f_moist + f_i_layer(i) + endif + enddo + + ! fraction of the atmosphere that is dry + f_dry = max(1.0_dp - f_moist, 1.0e-40_dp) + ! mixing ratios of dry species + do i = 1,d%sp%ng + if (d%sp_type(i) == DrySpeciesType) then + f_i_layer(i) = f_dry*d%f_i_dry(i) + endif + enddo + + end subroutine + + function dT_dP(d, P) result(dTdP) + type(AdiabatRCProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp) :: dTdP + + real(dp) :: P2, P1, deltaP, log10P + real(dp) :: T2, T1 + + log10P = log10(P) + P2 = log10P + log10P*d%epsj + P2 = max(min(P2, log10(d%P_surf)),log10(d%P_top)) + P2 = min(P2, log10(d%P_root)) + P1 = log10P - log10P*d%epsj + P1 = max(min(P1, log10(d%P_surf)),log10(d%P_top)) + deltaP = P2 - P1 + + call d%T%evaluate(P2, T2) + call d%T%evaluate(P1, T1) + + dTdP = (T2 - T1)/deltaP + + end function + + function general_adiabat_lapse_rate(d, P) result(dlnT_dlnP) + use clima_eqns, only: heat_capacity_eval + use clima_eqns_water, only: Rgas_cgs => Rgas + type(AdiabatRCProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp) :: dlnT_dlnP + + real(dp) :: T + real(dp) :: f_dry, cp_dry + real(dp) :: cp, L + real(dp) :: first_sumation, second_sumation, Beta_i + logical :: found + integer :: i + + real(dp), parameter :: Rgas_si = Rgas_cgs/1.0e7_dp ! ideal gas constant in SI units (J/(mol*K)) + + call d%T%evaluate(log10(P),T) + call mixing_ratios(d, P, T, d%f_i_cur, f_dry) + + ! heat capacity and latent heat + cp_dry = tiny(0.0_dp) + do i = 1,d%sp%ng + if (d%sp_type(i) == CondensingSpeciesType) then + L = d%sp%g(i)%sat%latent_heat(T) ! erg/g + L = L*d%sp%g(i)%mass*1.0e-7_dp ! convert to J/mol + d%L_i_cur(i) = L + endif + + ! heat capacity in J/(mol*K) + call heat_capacity_eval(d%sp%g(i)%thermo, T, found, cp) + if (.not. found) then + d%err = "Failed to compute heat capacity" + return + endif + d%cp_i_cur(i) = cp + if (d%sp_type(i) == DrySpeciesType) then + cp_dry = cp_dry + d%f_i_dry(i)*cp ! J/(mol*K) + endif + enddo + + ! Two sums in Equation 1 of Graham et al. (2021) + first_sumation = 0.0_dp + second_sumation = 0.0_dp + do i = 1,d%sp%ng + if (d%sp_type(i) == CondensingSpeciesType) then + Beta_i = d%L_i_cur(i)/(Rgas_si*T) + first_sumation = first_sumation + & + d%f_i_cur(i)*(d%cp_i_cur(i) - Rgas_si*Beta_i + Rgas_si*Beta_i**2.0_dp) + + second_sumation = second_sumation + & + Beta_i*d%f_i_cur(i) + endif + enddo + + ! Equation 1 in Graham et al. (2021), except simplified to assume no condensate present. + ! Units are SI + dlnT_dlnP = 1.0_dp/(f_dry*((cp_dry*f_dry + first_sumation)/(Rgas_si*(f_dry + second_sumation))) + second_sumation) + + end function + + subroutine right_hand_side(d, P, u, du) + use clima_eqns, only: heat_capacity_eval, gravity + use clima_eqns_water, only: Rgas_cgs => Rgas + type(AdiabatRCProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp), intent(in) :: u(:) + real(dp), intent(out) :: du(:) + + real(dp) :: z + real(dp) :: T, f_dry, mubar, grav, dz_dP + integer :: i + + ! unpack u + z = u(1) + + call d%T%evaluate(log10(P),T) + call mixing_ratios(d, P, T, d%f_i_cur, f_dry) + + mubar = 0.0_dp + do i = 1,d%sp%ng + mubar = mubar + d%f_i_cur(i)*d%sp%g(i)%mass + enddo + + grav = gravity(d%planet_radius, d%planet_mass, z) + dz_dP = -(Rgas_cgs*T)/(grav*P*mubar) + + du(:) = [dz_dP] + + end subroutine + +end module \ No newline at end of file