From b4ca5890d35d716cc42d8f74c7331498d6a5faad Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Wed, 21 Aug 2024 15:33:10 -0700 Subject: [PATCH] added dry routines to make atmosphere --- clima/cython/AdiabatClimate.pyx | 62 ++++ clima/cython/AdiabatClimate_pxd.pxd | 17 + clima/fortran/AdiabatClimate.f90 | 49 +++ src/CMakeLists.txt | 1 + src/adiabat/clima_adiabat.f90 | 98 ++++++ src/adiabat/clima_adiabat_dry.f90 | 352 +++++++++++++++++++ templates/AdiabatClimate/Earth/settings.yaml | 4 +- tests/test_adiabat.f90 | 27 +- 8 files changed, 604 insertions(+), 6 deletions(-) create mode 100644 src/adiabat/clima_adiabat_dry.f90 diff --git a/clima/cython/AdiabatClimate.pyx b/clima/cython/AdiabatClimate.pyx index 3c78d26..2fc3899 100644 --- a/clima/cython/AdiabatClimate.pyx +++ b/clima/cython/AdiabatClimate.pyx @@ -87,6 +87,37 @@ cdef class AdiabatClimate: &ng, P_i_surf.data, err) if len(err.strip()) > 0: raise ClimaException(err.decode("utf-8").strip()) + + def make_profile_dry(self, ndarray[double, ndim=1] P, ndarray[double, ndim=1] T, ndarray[double, ndim=2] f_i): + """Given a P, T and mixing ratios, this function will update all atmospheric variables + (except self%P_trop and self%convecting_with_below) to reflect these inputs. + The atmosphere is assumed to be dry (no condensibles). Any gas exceeding saturation + will not be altered. + + Parameters + ---------- + P : ndarray[double,ndim=1] + Pressure (dynes/cm^2). The first element is the surface. + T : ndarray[double,ndim=1] + Temperature (K) defined on `P`. + f_i : ndarray[double,ndim=2] + Mixing ratios defined on `P` of shape (size(P),ng). + """ + cdef int dim_P = P.shape[0] + cdef int dim_T = T.shape[0] + cdef ndarray f_i_ = np.asfortranarray(f_i) + cdef int dim1_f_i = f_i_.shape[0] + cdef int dim2_f_i = f_i_.shape[1] + cdef char err[ERR_LEN+1] + wa_pxd.adiabatclimate_make_profile_dry_wrapper( + self._ptr, + &dim_P, P.data, + &dim_T, T.data, + &dim1_f_i, &dim2_f_i, f_i_.data, + err + ) + if len(err.strip()) > 0: + raise ClimaException(err.decode("utf-8").strip()) def make_column(self, double T_surf, ndarray[double, ndim=1] N_i_surf): """Similar to `make_profile`, but instead the input is column reservoirs @@ -214,6 +245,37 @@ cdef class AdiabatClimate: if len(err.strip()) > 0: raise ClimaException(err.decode("utf-8").strip()) return ISR, OLR + + def TOA_fluxes_dry(self, ndarray[double, ndim=1] P, ndarray[double, ndim=1] T, ndarray[double, ndim=2] f_i): + """Calls `make_profile_dry`, then does radiative transfer on the constructed atmosphere + + Parameters + ---------- + P : ndarray[double,ndim=1] + Pressure (dynes/cm^2). The first element is the surface. + T : ndarray[double,ndim=1] + Temperature (K) defined on `P`. + f_i : ndarray[double,ndim=2] + Mixing ratios defined on `P` of shape (size(P),ng). + """ + cdef int dim_P = P.shape[0] + cdef int dim_T = T.shape[0] + cdef ndarray f_i_ = np.asfortranarray(f_i) + cdef int dim1_f_i = f_i_.shape[0] + cdef int dim2_f_i = f_i_.shape[1] + cdef char err[ERR_LEN+1] + cdef double ISR, OLR; + wa_pxd.adiabatclimate_toa_fluxes_dry_wrapper( + self._ptr, + &dim_P, P.data, + &dim_T, T.data, + &dim1_f_i, &dim2_f_i, f_i_.data, + &ISR, &OLR, + err + ) + if len(err.strip()) > 0: + raise ClimaException(err.decode("utf-8").strip()) + return ISR, OLR def surface_temperature(self, ndarray[double, ndim=1] P_i_surf, double T_guess = 280): """Does a non-linear solve for the surface temperature that balances incoming solar diff --git a/clima/cython/AdiabatClimate_pxd.pxd b/clima/cython/AdiabatClimate_pxd.pxd index 0834ce5..6e4e1a1 100644 --- a/clima/cython/AdiabatClimate_pxd.pxd +++ b/clima/cython/AdiabatClimate_pxd.pxd @@ -28,6 +28,14 @@ cdef extern void adiabatclimate_make_column_wrapper(AdiabatClimate *ptr, double cdef extern void adiabatclimate_make_profile_bg_gas_wrapper(AdiabatClimate *ptr, double *T_surf, int *ng, double *P_i_surf, double *P_surf, char *bg_gas, char *err) +cdef extern void adiabatclimate_make_profile_dry_wrapper( + AdiabatClimate *ptr, + int *dim_P, double *P, + int *dim_T, double *T, + int *dim1_f_i, int *dim2_f_i, double *f_i, + char *err +) + cdef extern void adiabatclimate_toa_fluxes_wrapper(AdiabatClimate *ptr, double *T_surf, int *ng, double *P_i_surf, double *ISR, double *OLR, char *err) @@ -38,6 +46,15 @@ cdef extern void adiabatclimate_toa_fluxes_bg_gas_wrapper(AdiabatClimate *ptr, d double *P_i_surf, double *P_surf, char *bg_gas, double *ISR, double *OLR, char *err) +cdef extern void adiabatclimate_toa_fluxes_dry_wrapper( + AdiabatClimate *ptr, + int *dim_P, double *P, + int *dim_T, double *T, + int *dim1_f_i, int *dim2_f_i, double *f_i, + double *ISR, double *OLR, + char *err +) + cdef extern void adiabatclimate_surface_temperature_wrapper(AdiabatClimate *ptr, int *ng, double *P_i_surf, double *T_guess, double *T_surf, char *err) diff --git a/clima/fortran/AdiabatClimate.f90 b/clima/fortran/AdiabatClimate.f90 index f9fb48b..d43b375 100644 --- a/clima/fortran/AdiabatClimate.f90 +++ b/clima/fortran/AdiabatClimate.f90 @@ -142,6 +142,30 @@ subroutine adiabatclimate_make_profile_bg_gas_wrapper(ptr, T_surf, ng, P_i_surf, end subroutine +subroutine adiabatclimate_make_profile_dry_wrapper(ptr, dim_P, P, dim_T, T, dim1_f_i, dim2_f_i, f_i, err) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), value, intent(in) :: ptr + integer(c_int), intent(in) :: dim_P + real(c_double), intent(in) :: P(dim_P) + integer(c_int), intent(in) :: dim_T + real(c_double), intent(in) :: T(dim_T) + integer(c_int), intent(in) :: dim1_f_i, dim2_f_i + real(c_double), intent(in) :: f_i(dim1_f_i,dim2_f_i) + character(c_char), intent(out) :: err(err_len+1) + + character(:), allocatable :: err_f + type(AdiabatClimate), pointer :: c + + call c_f_pointer(ptr, c) + call c%make_profile_dry(P, T, f_i, err_f) + + err(1) = c_null_char + if (allocated(err_f)) then + call copy_string_ftoc(err_f, err) + endif + +end subroutine + subroutine adiabatclimate_toa_fluxes_wrapper(ptr, T_surf, ng, P_i_surf, ISR, OLR, err) bind(c) use clima, only: AdiabatClimate type(c_ptr), value, intent(in) :: ptr @@ -217,6 +241,31 @@ subroutine adiabatclimate_toa_fluxes_bg_gas_wrapper(ptr, T_surf, ng, P_i_surf, P end subroutine +subroutine adiabatclimate_toa_fluxes_dry_wrapper(ptr, dim_P, P, dim_T, T, dim1_f_i, dim2_f_i, f_i, ISR, OLR, err) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), value, intent(in) :: ptr + integer(c_int), intent(in) :: dim_P + real(c_double), intent(in) :: P(dim_P) + integer(c_int), intent(in) :: dim_T + real(c_double), intent(in) :: T(dim_T) + integer(c_int), intent(in) :: dim1_f_i, dim2_f_i + real(c_double), intent(in) :: f_i(dim1_f_i,dim2_f_i) + real(c_double), intent(out) :: ISR, OLR + character(c_char), intent(out) :: err(err_len+1) + + character(:), allocatable :: err_f + type(AdiabatClimate), pointer :: c + + call c_f_pointer(ptr, c) + call c%TOA_fluxes_dry(P, T, f_i, ISR, OLR, err_f) + + err(1) = c_null_char + if (allocated(err_f)) then + call copy_string_ftoc(err_f, err) + endif + +end subroutine + subroutine adiabatclimate_surface_temperature_wrapper(ptr, ng, P_i_surf, T_guess, T_surf, err) bind(c) use clima, only: AdiabatClimate type(c_ptr), value, intent(in) :: ptr diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 458b34e..cd01b2d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,6 +33,7 @@ set(ADIABAT_SOURCES adiabat/clima_adiabat_general.f90 adiabat/clima_adiabat_rc.f90 adiabat/clima_adiabat_solve.f90 + adiabat/clima_adiabat_dry.f90 adiabat/clima_adiabat.f90 ) diff --git a/src/adiabat/clima_adiabat.f90 b/src/adiabat/clima_adiabat.f90 index 4996047..316d867 100644 --- a/src/adiabat/clima_adiabat.f90 +++ b/src/adiabat/clima_adiabat.f90 @@ -130,10 +130,12 @@ module clima_adiabat procedure :: make_profile => AdiabatClimate_make_profile procedure :: make_column => AdiabatClimate_make_column procedure :: make_profile_bg_gas => AdiabatClimate_make_profile_bg_gas + procedure :: make_profile_dry => AdiabatClimate_make_profile_dry ! Constructs atmosphere and does radiative transfer procedure :: TOA_fluxes => AdiabatClimate_TOA_fluxes procedure :: TOA_fluxes_column => AdiabatClimate_TOA_fluxes_column procedure :: TOA_fluxes_bg_gas => AdiabatClimate_TOA_fluxes_bg_gas + procedure :: TOA_fluxes_dry => AdiabatClimate_TOA_fluxes_dry ! Non-linear solves for equilibrium climate procedure :: surface_temperature => AdiabatClimate_surface_temperature procedure :: surface_temperature_column => AdiabatClimate_surface_temperature_column @@ -493,6 +495,76 @@ subroutine fcn(n_, x_, fvec_, iflag_) end subroutine end subroutine + !> Given a P, T and mixing ratios, this function will update all atmospheric variables + !> (except self%P_trop and self%convecting_with_below) to reflect these inputs. + !> The atmosphere is assumed to be dry (no condensibles). Any gas exceeding saturation + !> will not be altered. + subroutine AdiabatClimate_make_profile_dry(self, P, T, f_i, err) + use clima_adiabat_dry, only: make_profile_dry + use clima_const, only: k_boltz, N_avo + class(AdiabatClimate), intent(inout) :: self + real(dp), intent(in) :: P(:) !! Pressure (dynes/cm^2). The first element is the surface. + real(dp), intent(in) :: T(:) !! Temperature (K) defined on `P`. + real(dp), intent(in) :: f_i(:,:) !! Mixing ratios defined on `P` of shape (size(P),ng). + character(:), allocatable, intent(out) :: err + + real(dp), allocatable :: P_e(:), z_e(:), T_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),T_e(2*self%nz+1)) + allocate(f_i_e(2*self%nz+1,self%sp%ng),lapse_rate_e(2*self%nz+1)) + allocate(density(self%nz)) + + call make_profile_dry(P, T, f_i, & + self%sp, self%nz, & + self%planet_mass, self%planet_radius, self%P_top, & + self%rtol, self%atol, & + P_e, z_e, T_e, f_i_e, lapse_rate_e, & + err) + if (allocated(err)) return + + ! We assume the atmosphere is dry + self%N_surface = 0.0_dp + self%N_ocean = 0.0_dp + + self%T_surf = T_e(1) + self%P_surf = P_e(1) + do i = 1,self%nz + self%P(i) = P_e(2*i) + self%T(i) = T_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 + + ! lapse rates + self%lapse_rate_intended(1) = lapse_rate_e(1) + do i = 2,self%nz + self%lapse_rate_intended(i) = lapse_rate_e(2*i-2) + enddo + + ! Convecting with below is not changed. We don't get convecting information. + + self%lapse_rate(1) = (log(self%T(1)) - log(self%T_surf))/(log(self%P(1)) - log(self%P_surf)) + do i = 2,self%nz + self%lapse_rate(i) = (log(self%T(i)) - log(self%T(i-1)))/(log(self%P(i)) - log(self%P(i-1))) + enddo + + end subroutine + !> Copies the atmosphere to the data used for radiative transfer subroutine AdiabatClimate_copy_atm_to_radiative_grid(self) class(AdiabatClimate), intent(inout) :: self @@ -597,6 +669,32 @@ subroutine AdiabatClimate_TOA_fluxes_bg_gas(self, T_surf, P_i_surf, P_surf, bg_g if (allocated(err)) return end subroutine + + !> Calls `make_profile_dry`, then does radiative transfer on the constructed atmosphere + subroutine AdiabatClimate_TOA_fluxes_dry(self, P, T, f_i, ISR, OLR, err) + class(AdiabatClimate), intent(inout) :: self + real(dp), intent(in) :: P(:) !! Pressure (dynes/cm^2). The first element is the surface. + real(dp), intent(in) :: T(:) !! Temperature (K) defined on `P`. + real(dp), intent(in) :: f_i(:,:) !! Mixing ratios defined on `P` of shape (size(P),ng). + real(dp), intent(out) :: ISR !! Top-of-atmosphere incoming solar radiation (mW/m^2) + real(dp), intent(out) :: OLR !! Top-of-atmosphere outgoing longwave radiation (mW/m^2) + character(:), allocatable, intent(out) :: err + + ! make atmosphere profile + call self%make_profile_dry(P, T, f_i, err) + if (allocated(err)) return + + ! set albedo + if (associated(self%albedo_fcn)) then + self%rad%surface_albedo = self%albedo_fcn(self%T_surf) + endif + + ! Do radiative transfer + call self%copy_atm_to_radiative_grid() + call self%rad%TOA_fluxes(self%T_surf, self%T_r, self%P_r/1.0e6_dp, self%densities_r, self%dz_r, ISR=ISR, OLR=OLR, err=err) + if (allocated(err)) return + + end subroutine !> Does a non-linear solve for the surface temperature that balances incoming solar !> and outgoing longwave radiation. Uses `make_profile`. diff --git a/src/adiabat/clima_adiabat_dry.f90 b/src/adiabat/clima_adiabat_dry.f90 new file mode 100644 index 0000000..64ad501 --- /dev/null +++ b/src/adiabat/clima_adiabat_dry.f90 @@ -0,0 +1,352 @@ +module clima_adiabat_dry + 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 :: AdiabatDryProfileData + ! Input data + real(dp), pointer :: P_in(:) + real(dp), pointer :: T_in(:) + real(dp), pointer :: f_i_in(:,:) + 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 :: rtol + real(dp), pointer :: atol + ! Ouput + real(dp), pointer :: P(:) + real(dp), pointer :: z(:) + real(dp), pointer :: T(:) + real(dp), pointer :: f_i(:,:) + real(dp), pointer :: lapse_rate(:) + + real(dp) :: P_surf !! surface pressure from inputs + + type(linear_interp_1d) :: T_interp !! Interpolates input T + type(linear_interp_1d), allocatable :: f_i_interp(:) !! Interpolates mixing ratios (ng) + + !> index for saving T, z and mixing ratios + integer :: j + + !> work variables. All dimension (ng) + real(dp), allocatable :: f_i_cur(:) + + !> This helps us propogate error messages outside of the integration + character(:), allocatable :: err + + end type + + type, extends(dop853_class) :: dop853_custom + type(AdiabatDryProfileData), pointer :: d => NULL() + end type + +contains + + subroutine make_profile_dry(P_in, T_in, f_i_in, & + sp, nz, & + planet_mass, planet_radius, P_top, & + rtol, atol, & + P, z, T, f_i, lapse_rate, & + err) + use futils, only: linspace + real(dp), target, intent(in) :: P_in(:) + real(dp), target, intent(in) :: T_in(:) + real(dp), target, intent(in) :: f_i_in(:,:) + 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 + real(dp), target, intent(in) :: rtol, atol + real(dp), target, intent(out) :: P(:), z(:), T(:) + real(dp), target, intent(out) :: f_i(:,:) + real(dp), target, intent(out) :: lapse_rate(:) + character(:), allocatable, intent(out) :: err + + real(dp), allocatable :: tmp_arr(:) + real(dp), allocatable :: log10P_in_interp(:) + type(AdiabatDryProfileData) :: d + integer :: i, ierr + + ! Check inputs + if (any(T_in < 0.0_dp)) then + err = '`T` can not have negative elements' + return + endif + if (any(P_in < 0.0_dp)) then + err = '`P` can not have negative elements' + return + endif + if (P_in(1) < P_top) then + err = 'The first element of `P` must be greater than `P_top`' + return + endif + if (size(T_in) <= 1) then + err = '`P` must have a size > 1.' + return + endif + if (size(T_in) /= size(P_in)) then + err = '`T` and `P` must have the same length' + return + endif + do i = 1,size(P_in)-1 + if (P_in(i+1) >= P_in(i)) then + err = '`P` must be strictly decreasing' + return + endif + enddo + if (any(f_i_in < 0.0_dp)) then + err = '`f_i` can not have negative elements' + return + endif + if (size(f_i_in,1) /= size(P_in)) then + err = 'The first dimension of `f_i` must match the size of `P`' + return + endif + if (size(f_i_in,2) /= sp%ng) then + err = 'The second dimension of `f_i` must match the the number of species' + return + endif + + ! Associate inputs + d%P_in => P_in + d%T_in => T_in + d%f_i_in => f_i_in + d%sp => sp + d%nz => nz + d%planet_mass => planet_mass + d%planet_radius => planet_radius + d%P_top => P_top + d%rtol => rtol + d%atol => atol + ! outputs + d%P => P + d%z => z + d%T => T + d%f_i => f_i + d%lapse_rate => lapse_rate + + ! allocate + allocate(d%f_i_cur(sp%ng)) + + ! Make P profile + d%P_surf = P_in(1) + 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 + + ! Build the interpolators + allocate(log10P_in_interp(size(P_in))) + log10P_in_interp = log10(P_in) + log10P_in_interp = log10P_in_interp(size(log10P_in_interp):1:-1) + + ! Temperature + allocate(tmp_arr(size(T_in))) + tmp_arr = T_in + tmp_arr = tmp_arr(size(tmp_arr):1:-1) + call d%T_interp%initialize(log10P_in_interp, tmp_arr, ierr) + if (ierr /= 0) then + err = 'Failed to initialize interpolator in "make_profile_dry"' + return + endif + + ! Mixing ratios + allocate(d%f_i_interp(sp%ng)) + do i = 1,sp%ng + tmp_arr = log10(f_i_in(:,i)) + tmp_arr = tmp_arr(size(tmp_arr):1:-1) + call d%f_i_interp(i)%initialize(log10P_in_interp, tmp_arr, ierr) + if (ierr /= 0) then + err = 'Failed to initialize interpolator in "make_profile_dry"' + return + endif + enddo + + call integrate(d, err) + if (allocated(err)) return + + end subroutine + + subroutine integrate(d, err) + type(AdiabatDryProfileData), target, intent(inout) :: d + character(:), allocatable, intent(out) :: err + + type(dop853_custom) :: dop + logical :: status_ok + integer :: idid + real(dp) :: Pn, u(1) + + 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 + call mixing_ratios(d, d%P_surf, d%f_i(1,:)) + call d%T_interp%evaluate(log10(d%P_surf),d%T(1)) + d%z(1) = 0.0_dp + call dry_adiabat(d, d%P_surf, d%lapse_rate(1), err) + if (allocated(err)) return + if (.not. status_ok) then + err = 'dop853 initialization failed' + return + endif + + u = [0.0_dp] + Pn = d%P_surf + 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 + + 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(AdiabatDryProfileData), pointer :: d + real(dp) :: z_cur, P_cur, P_old + + P_old = xold + P_cur = x + z_cur = y(1) + + select type (self) + class is (dop853_custom) + 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 + + ! save the results + if (d%j <= size(d%P)) then + if (d%P(d%j) <= P_old .and. d%P(d%j) >= P_cur) then + do while (d%P(d%j) >= P_cur) + d%z(d%j) = self%contd8(1, d%P(d%j)) + call d%T_interp%evaluate(log10(d%P(d%j)), d%T(d%j)) + call mixing_ratios(d, d%P(d%j), d%f_i(d%j,:)) + call dry_adiabat(d, d%P(d%j), d%lapse_rate(d%j), d%err) + if (allocated(d%err)) return + d%j = d%j + 1 + if (d%j > size(d%P)) exit + enddo + endif + endif + + 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_custom) + call right_hand_side(self%d, P, u, du) + end select + + end subroutine + + subroutine dry_adiabat(d, P, lapse_rate, err) + use clima_eqns, only: heat_capacity_eval + use clima_eqns_water, only: Rgas_cgs => Rgas + type(AdiabatDryProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp), intent(out) :: lapse_rate + character(:), allocatable, intent(out) :: err + + real(dp) :: T, cp_tmp, cp + logical :: found + integer :: i + + real(dp), parameter :: Rgas_si = Rgas_cgs/1.0e7_dp ! ideal gas constant in SI units (J/(mol*K)) + + ! Get temperature + call d%T_interp%evaluate(log10(P), T) + + ! Get mixing ratios + call mixing_ratios(d, P, d%f_i_cur) + + cp = tiny(0.0_dp) + do i = 1,d%sp%ng + call heat_capacity_eval(d%sp%g(i)%thermo, T, 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*d%f_i_cur(i) ! J/(mol*K) + enddo + + lapse_rate = Rgas_si/cp + + end subroutine + + subroutine mixing_ratios(d, P, f_i_layer) + type(AdiabatDryProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp), intent(out) :: f_i_layer(:) + + integer :: i + + do i = 1,d%sp%ng + call d%f_i_interp(i)%evaluate(log10(P),f_i_layer(i)) + f_i_layer(i) = 10.0_dp**f_i_layer(i) + enddo + + end subroutine + + subroutine right_hand_side(d, P, u, du) + use clima_eqns, only: gravity + use clima_eqns_water, only: Rgas_cgs => Rgas + type(AdiabatDryProfileData), intent(inout) :: d + real(dp), intent(in) :: P + real(dp), intent(in) :: u(:) + real(dp), intent(out) :: du(:) + + real(dp) :: z + real(dp) :: T, mubar, grav, dz_dP + integer :: i + + ! unpack u + z = u(1) + + ! Get temperature + call d%T_interp%evaluate(log10(P), T) + + ! Get mixing ratios + call mixing_ratios(d, P, d%f_i_cur) + + mubar = 0.0_dp + do i = 1,d%sp%ng + mubar = mubar + d%f_i_cur(i)*d%sp%g(i)%mass + enddo + + ! rate of change of altitude + 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 diff --git a/templates/AdiabatClimate/Earth/settings.yaml b/templates/AdiabatClimate/Earth/settings.yaml index a6b9767..a3c4c00 100644 --- a/templates/AdiabatClimate/Earth/settings.yaml +++ b/templates/AdiabatClimate/Earth/settings.yaml @@ -28,7 +28,7 @@ optical-properties: # If you choose method RandomOverlapResortRebin, then you must supply # `number-of-bins`, which is the number of bins the method rebins to - number-of-bins: 32 + number-of-bins: 16 # Here you specify opacites. You can see avaliable opacities by looking in # `clima/data/` folder. `CIA`, `rayleigh`, `photolysis-xs` and `water-continuum` @@ -39,5 +39,5 @@ optical-properties: opacities: {k-distributions: on, CIA: on, rayleigh: on} solar: k-method: RandomOverlapResortRebin - number-of-bins: 32 + number-of-bins: 16 opacities: {k-distributions: on, CIA: on, rayleigh: on} diff --git a/tests/test_adiabat.f90 b/tests/test_adiabat.f90 index c432e0e..8cf1182 100644 --- a/tests/test_adiabat.f90 +++ b/tests/test_adiabat.f90 @@ -6,8 +6,8 @@ program test type(AdiabatClimate) :: c character(:), allocatable :: err - real(dp) :: T, OLR - real(dp), allocatable :: P_i_surf(:), N_i_surf(:) + real(dp) :: T, OLR, ISR, OLR1, ISR1 + real(dp), allocatable :: P_i_surf(:), N_i_surf(:), f_i(:,:) integer :: i logical :: converged procedure(ocean_solubility_fcn), pointer :: ocean_fcn_ptr @@ -71,14 +71,33 @@ program test c%T_trop = 150.0_dp P_i_surf = [270.0_dp, 10.0_dp, 1.0_dp, 1.0e-10_dp, 1.0e-10_dp, 1.0e-10_dp, 1.0e-10_dp] P_i_surf = P_i_surf*1.0e6_dp - call c%make_profile(280.0_dp, & - P_i_surf, & + call c%TOA_fluxes(280.0_dp, & + P_i_surf, ISR, OLR, & err=err) if (allocated(err)) then print*,err stop 1 endif + ! Test make_profile_dry + allocate(f_i(c%nz+1,c%sp%ng)) + f_i(1,:) = c%f_i(1,:) + do i = 1,c%nz + f_i(i+1,:) = c%f_i(i,:) + enddo + call c%TOA_fluxes_dry( & + [c%P_surf,c%P], & + [c%T_surf,c%T], & + f_i, & + ISR1, OLR1, & + err & + ) + if (allocated(err)) then + print*,err + stop 1 + endif + print*,ISR1/ISR, OLR1/OLR + ! Test ocean solubility functionality ocean_fcn_ptr => ocean_fcn call c%set_ocean_solubility_fcn('H2O', ocean_fcn_ptr, err)