From ab8e44f322adc5bdb651afde84fc7fabfff0652f Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Wed, 21 Feb 2024 15:19:36 -0800 Subject: [PATCH] python wrapper. clean up --- clima/cython/AdiabatClimate.pyx | 147 +++++++++++++++++-- clima/cython/AdiabatClimate_pxd.pxd | 33 ++++- clima/cython/_clima.pyx | 4 +- clima/fortran/AdiabatClimate.f90 | 211 +++++++++++++++++++++++++++- src/adiabat/clima_adiabat.f90 | 19 ++- src/adiabat/clima_adiabat_rc.f90 | 10 ++ src/adiabat/clima_adiabat_solve.f90 | 100 +++++++++---- 7 files changed, 475 insertions(+), 49 deletions(-) diff --git a/clima/cython/AdiabatClimate.pyx b/clima/cython/AdiabatClimate.pyx index 761e232..4696dd6 100644 --- a/clima/cython/AdiabatClimate.pyx +++ b/clima/cython/AdiabatClimate.pyx @@ -12,7 +12,7 @@ cdef class AdiabatClimate: cdef void *_ptr def __init__(self, species_file = None, settings_file = None, - flux_file = None, data_dir = None): + flux_file = None, data_dir = None, cbool double_radiative_grid = True): """Initializes `AdiabatClimate` Parameters @@ -47,8 +47,8 @@ cdef class AdiabatClimate: # Initialize wa_pxd.adiabatclimate_create_wrapper(&self._ptr, species_file_c, - settings_file_c, flux_file_c, data_dir_c, - err) + settings_file_c, flux_file_c, data_dir_c, &double_radiative_grid, + err) if len(err.strip()) > 0: raise ClimaException(err.decode("utf-8").strip()) @@ -316,6 +316,47 @@ cdef class AdiabatClimate: raise ClimaException(err.decode("utf-8").strip()) return T_surf + def RCE(self, ndarray[double, ndim=1] P_i_surf, double T_surf_guess, ndarray[double, ndim=1] T_guess, + convecting_with_below = None): + """Compute full radiative-convective equilibrium. + + Parameters + ---------- + P_i_surf : ndarray[double,ndim=1] + Array of surface pressures of each species (dynes/cm^2) + T_surf_guess : float + A guess for the surface temperature (K) + T_guess : ndarray[double,ndim=1] + A guess for the temperature in each atmospheric layer (K) + convecting_with_below : ndarray[bool,ndim=1], optional + An array describing a guess for the radiative vs. convective + regions of the atmosphere + + Returns + ------- + bool + Whether the routine converged or not. + """ + cdef int ng = P_i_surf.shape[0] + cdef int dim_T_guess = T_guess.shape[0] + cdef ndarray[cbool, ndim=1] convecting_with_below_ = np.array([False],dtype=np.bool_) + cdef cbool convecting_with_below_present + if convecting_with_below is None: + convecting_with_below_present = False + else: + convecting_with_below_present = True + convecting_with_below_ = convecting_with_below + cdef int dim_convecting_with_below = convecting_with_below_.shape[0] + cdef cbool converged + cdef char err[ERR_LEN+1] + + wa_pxd.adiabatclimate_rce_wrapper(&self._ptr, &ng, P_i_surf.data, &T_surf_guess, + &dim_T_guess, T_guess.data, &convecting_with_below_present, + &dim_convecting_with_below, convecting_with_below_.data, &converged, err) + if len(err.strip()) > 0: + raise ClimaException(err.decode("utf-8").strip()) + return converged + def to_regular_grid(self): "Re-grids atmosphere so that each grid cell is equally spaced in altitude." cdef char err[ERR_LEN+1] @@ -323,7 +364,7 @@ cdef class AdiabatClimate: if len(err.strip()) > 0: raise ClimaException(err.decode("utf-8").strip()) - def out2atmosphere_txt(self, filename, ndarray[double, ndim=1] eddy, bool overwrite = False, bool clip = True): + def out2atmosphere_txt(self, filename, ndarray[double, ndim=1] eddy, cbool overwrite = False, cbool clip = True): """Saves state of the atmosphere to a file. Parameters @@ -390,10 +431,10 @@ cdef class AdiabatClimate: use the initial guess in `self.make_column_P_guess` """ def __get__(self): - cdef bool val + cdef cbool val wa_pxd.adiabatclimate_use_make_column_p_guess_get(&self._ptr, &val) return val - def __set__(self, bool val): + def __set__(self, cbool val): wa_pxd.adiabatclimate_use_make_column_p_guess_set(&self._ptr, &val) property make_column_P_guess: @@ -418,10 +459,10 @@ cdef class AdiabatClimate: it matches the skin temperature. The initial guess will always be self.T_trop. """ def __get__(self): - cdef bool val + cdef cbool val wa_pxd.adiabatclimate_solve_for_t_trop_get(&self._ptr, &val) return val - def __set__(self, bool val): + def __set__(self, cbool val): wa_pxd.adiabatclimate_solve_for_t_trop_set(&self._ptr, &val) property albedo_fcn: @@ -478,10 +519,10 @@ cdef class AdiabatClimate: observed dayside temperature of a tidally locked planet. """ def __get__(self): - cdef bool val + cdef cbool val wa_pxd.adiabatclimate_tidally_locked_dayside_get(&self._ptr, &val) return val - def __set__(self, bool val): + def __set__(self, cbool val): wa_pxd.adiabatclimate_tidally_locked_dayside_set(&self._ptr, &val) property L: @@ -546,7 +587,46 @@ cdef class AdiabatClimate: var = Radtran() var._ptr = ptr1 return var + + property convecting_with_below: + """ndarray[bool,ndim=1], shape (nz). If True, then the layer below + is convecting with the current layer. Index 1 determines if the + first atomspheric layer is convecting with the ground. + """ + def __get__(self): + cdef int dim1 + wa_pxd.adiabatclimate_convecting_with_below_get_size(&self._ptr, &dim1) + cdef ndarray[cbool, ndim=1] arr = np.empty(dim1, bool) + wa_pxd.adiabatclimate_convecting_with_below_get(&self._ptr, &dim1, arr.data) + return arr + + property lapse_rate: + "ndarray[double,ndim=1], shape (nz). The true lapse rate (dlnT/dlnP)." + def __get__(self): + cdef int dim1 + wa_pxd.adiabatclimate_lapse_rate_get_size(&self._ptr, &dim1) + cdef ndarray arr = np.empty(dim1, np.double) + wa_pxd.adiabatclimate_lapse_rate_get(&self._ptr, &dim1, arr.data) + return arr + + property lapse_rate_intended: + "ndarray[double,ndim=1], shape (nz). The computed lapse rate (dlnT/dlnP)." + def __get__(self): + cdef int dim1 + wa_pxd.adiabatclimate_lapse_rate_intended_get_size(&self._ptr, &dim1) + cdef ndarray arr = np.empty(dim1, np.double) + wa_pxd.adiabatclimate_lapse_rate_intended_get(&self._ptr, &dim1, arr.data) + return arr + property convective_newton_step_size: + "float. The size of the newton step." + def __get__(self): + cdef double val + wa_pxd.adiabatclimate_convective_newton_step_size_get(&self._ptr, &val) + return val + def __set__(self, double val): + wa_pxd.adiabatclimate_convective_newton_step_size_set(&self._ptr, &val) + property rtol: "float. Relative tolerance of integration." def __get__(self): @@ -573,6 +653,53 @@ cdef class AdiabatClimate: return val def __set__(self, double val): wa_pxd.adiabatclimate_tol_make_column_set(&self._ptr, &val) + + property epsj: + "float. Perturbation for the jacobian" + def __get__(self): + cdef double val + wa_pxd.adiabatclimate_epsj_get(&self._ptr, &val) + return val + def __set__(self, double val): + wa_pxd.adiabatclimate_epsj_set(&self._ptr, &val) + + property xtol_rc: + "float. xtol for RC equilibrium" + def __get__(self): + cdef double val + wa_pxd.adiabatclimate_xtol_rc_get(&self._ptr, &val) + return val + def __set__(self, double val): + wa_pxd.adiabatclimate_xtol_rc_set(&self._ptr, &val) + + property max_rc_iters: + "int. Max number of iterations in the RCE routine" + def __get__(self): + cdef int val + wa_pxd.adiabatclimate_max_rc_iters_get(&self._ptr, &val) + return val + def __set__(self, int val): + wa_pxd.adiabatclimate_max_rc_iters_set(&self._ptr, &val) + + property max_rc_iters_convection: + """int. Max number of iterations for which convective layers can + be converged to radiative layers in the RCE routine + """ + def __get__(self): + cdef int val + wa_pxd.adiabatclimate_max_rc_iters_convection_get(&self._ptr, &val) + return val + def __set__(self, int val): + wa_pxd.adiabatclimate_max_rc_iters_convection_set(&self._ptr, &val) + + property verbose: + "bool. verbosity" + def __get__(self): + cdef cbool val + wa_pxd.adiabatclimate_verbose_get(&self._ptr, &val) + return val + def __set__(self, cbool val): + wa_pxd.adiabatclimate_verbose_set(&self._ptr, &val) property P_surf: "float. Surface pressure (dynes/cm^2)" diff --git a/clima/cython/AdiabatClimate_pxd.pxd b/clima/cython/AdiabatClimate_pxd.pxd index 3a30724..050f16f 100644 --- a/clima/cython/AdiabatClimate_pxd.pxd +++ b/clima/cython/AdiabatClimate_pxd.pxd @@ -12,7 +12,7 @@ cdef extern void deallocate_adiabatclimate(void *ptr); # wrappers for functions cdef extern void adiabatclimate_create_wrapper(void *ptr, char *species_file, - char *settings_file, char *flux_file, char *data_dir, char *err); + char *settings_file, char *flux_file, char *data_dir, bool *double_radiative_grid, char *err); cdef extern void adiabatclimate_make_profile_wrapper(void *ptr, double *T_surf, int *ng, double *P_i_surf, char *err) @@ -43,6 +43,10 @@ cdef extern void adiabatclimate_surface_temperature_bg_gas_wrapper(void *ptr, in double *P_i_surf, double *P_surf, char *bg_gas, double *T_guess, double *T_surf, char *err) +cdef extern void adiabatclimate_rce_wrapper(void *ptr, int *ng, double *P_i_surf, double *T_surf_guess, + int *dim_T_guess, double *T_guess, bool *convecting_with_below_present, + int *dim_convecting_with_below, bool *convecting_with_below, bool *converged, char *err) + cdef extern void adiabatclimate_set_ocean_solubility_fcn_wrapper(void *ptr, char *species_c, ocean_solubility_fcn fcn, char *err) cdef extern void adiabatclimate_to_regular_grid_wrapper(void *ptr, char *err) @@ -99,6 +103,18 @@ cdef extern void adiabatclimate_species_names_get(void *ptr, int *dim1, char* sp cdef extern void adiabatclimate_rad_get(void *ptr, void *ptr1) +cdef extern void adiabatclimate_convecting_with_below_get_size(void *ptr, int *dim1) +cdef extern void adiabatclimate_convecting_with_below_get(void *ptr, int *dim1, bool *arr) + +cdef extern void adiabatclimate_lapse_rate_get_size(void *ptr, int *dim1) +cdef extern void adiabatclimate_lapse_rate_get(void *ptr, int *dim1, double *arr) + +cdef extern void adiabatclimate_lapse_rate_intended_get_size(void *ptr, int *dim1) +cdef extern void adiabatclimate_lapse_rate_intended_get(void *ptr, int *dim1, double *arr) + +cdef extern void adiabatclimate_convective_newton_step_size_get(void *ptr, double *val) +cdef extern void adiabatclimate_convective_newton_step_size_set(void *ptr, double *val) + cdef extern void adiabatclimate_rtol_get(void *ptr, double *val) cdef extern void adiabatclimate_rtol_set(void *ptr, double *val) @@ -108,6 +124,21 @@ cdef extern void adiabatclimate_atol_set(void *ptr, double *val) cdef extern void adiabatclimate_tol_make_column_get(void *ptr, double *val) cdef extern void adiabatclimate_tol_make_column_set(void *ptr, double *val) +cdef extern void adiabatclimate_epsj_get(void *ptr, double *val) +cdef extern void adiabatclimate_epsj_set(void *ptr, double *val) + +cdef extern void adiabatclimate_xtol_rc_get(void *ptr, double *val) +cdef extern void adiabatclimate_xtol_rc_set(void *ptr, double *val) + +cdef extern void adiabatclimate_max_rc_iters_get(void *ptr, int *val) +cdef extern void adiabatclimate_max_rc_iters_set(void *ptr, int *val) + +cdef extern void adiabatclimate_max_rc_iters_convection_get(void *ptr, int *val) +cdef extern void adiabatclimate_max_rc_iters_convection_set(void *ptr, int *val) + +cdef extern void adiabatclimate_verbose_get(void *ptr, bool *val) +cdef extern void adiabatclimate_verbose_set(void *ptr, bool *val) + cdef extern void adiabatclimate_p_surf_get(void *ptr, double *val) cdef extern void adiabatclimate_p_trop_get(void *ptr, double *val) diff --git a/clima/cython/_clima.pyx b/clima/cython/_clima.pyx index 5b76ea9..522a3d4 100644 --- a/clima/cython/_clima.pyx +++ b/clima/cython/_clima.pyx @@ -1,6 +1,6 @@ -from numpy cimport ndarray +from numpy cimport ndarray, uint8_t from libc.stdint cimport uintptr_t -from libcpp cimport bool +from libcpp cimport bool as cbool import numpy as np import ctypes as ct import os diff --git a/clima/fortran/AdiabatClimate.f90 b/clima/fortran/AdiabatClimate.f90 index 949309c..74c4bb2 100644 --- a/clima/fortran/AdiabatClimate.f90 +++ b/clima/fortran/AdiabatClimate.f90 @@ -26,19 +26,21 @@ subroutine deallocate_adiabatclimate(ptr) bind(c) !!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine adiabatclimate_create_wrapper(ptr, species_file, & - settings_file, flux_file, data_dir, err) bind(c) + settings_file, flux_file, data_dir, double_radiative_grid, err) bind(c) use clima, only: AdiabatClimate type(c_ptr), intent(in) :: ptr character(kind=c_char), intent(in) :: species_file(*) character(kind=c_char), intent(in) :: settings_file(*) character(kind=c_char), intent(in) :: flux_file(*) character(kind=c_char), intent(in) :: data_dir(*) + logical(c_bool), optional, intent(in) :: double_radiative_grid character(kind=c_char), intent(out) :: err(err_len+1) character(len=:), allocatable :: species_file_f character(len=:), allocatable :: settings_file_f character(len=:), allocatable :: flux_file_f character(len=:), allocatable :: data_dir_f + logical :: double_radiative_grid_f character(:), allocatable :: err_f type(AdiabatClimate), pointer :: c @@ -48,6 +50,7 @@ subroutine adiabatclimate_create_wrapper(ptr, species_file, & allocate(character(len=len_cstring(settings_file))::settings_file_f) allocate(character(len=len_cstring(flux_file))::flux_file_f) allocate(character(len=len_cstring(data_dir))::data_dir_f) + double_radiative_grid_f = double_radiative_grid call copy_string_ctof(species_file, species_file_f) call copy_string_ctof(settings_file, settings_file_f) @@ -58,6 +61,7 @@ subroutine adiabatclimate_create_wrapper(ptr, species_file, & settings_file_f, & flux_file_f, & data_dir_f, & + double_radiative_grid_f, & err_f) err(1) = c_null_char @@ -288,6 +292,46 @@ subroutine adiabatclimate_surface_temperature_bg_gas_wrapper(ptr, ng, P_i_surf, end subroutine +subroutine adiabatclimate_rce_wrapper(ptr, ng, P_i_surf, T_surf_guess, dim_T_guess, T_guess, & + convecting_with_below_present, dim_convecting_with_below, & + convecting_with_below, converged, err) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(in) :: ng + real(c_double), intent(in) :: P_i_surf(ng) + real(c_double), intent(in) :: T_surf_guess + integer(c_int), intent(in) :: dim_T_guess + real(c_double), intent(in) :: T_guess(dim_T_guess) + logical(c_bool), intent(in) :: convecting_with_below_present + integer(c_int), intent(in) :: dim_convecting_with_below + logical(c_bool), intent(in) :: convecting_with_below(dim_convecting_with_below) + logical(c_bool), intent(out) :: converged + character(c_char), intent(out) :: err(err_len+1) + + logical, allocatable :: convecting_with_below_f(:) + logical :: converged_f + character(:), allocatable :: err_f + type(AdiabatClimate), pointer :: c + + call c_f_pointer(ptr, c) + + allocate(convecting_with_below_f(dim_convecting_with_below)) + convecting_with_below_f = convecting_with_below + + converged_f = .false. + if (convecting_with_below_present) then + converged_f = c%RCE(P_i_surf, T_surf_guess, T_guess, convecting_with_below_f, err_f) + else + converged_f = c%RCE(P_i_surf, T_surf_guess, T_guess, err=err_f) + endif + converged = converged_f + err(1) = c_null_char + if (allocated(err_f)) then + call copy_string_ftoc(err_f, err) + endif + +end subroutine + subroutine adiabatclimate_set_ocean_solubility_fcn_wrapper(ptr, species_c, fcn_c, err) bind(c) use clima, only: AdiabatClimate use clima, only: ocean_solubility_fcn @@ -696,6 +740,81 @@ subroutine adiabatclimate_rad_get(ptr, ptr1) bind(c) ptr1 = c_loc(c%rad) end subroutine +subroutine adiabatclimate_convecting_with_below_get_size(ptr, dim1) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(out) :: dim1 + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + dim1 = size(c%convecting_with_below,1) +end subroutine + +subroutine adiabatclimate_convecting_with_below_get(ptr, dim1, arr) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(in) :: dim1 + logical(c_bool), intent(out) :: arr(dim1) + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + arr = c%convecting_with_below +end subroutine + +subroutine adiabatclimate_lapse_rate_get_size(ptr, dim1) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(out) :: dim1 + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + dim1 = size(c%lapse_rate,1) +end subroutine + +subroutine adiabatclimate_lapse_rate_get(ptr, dim1, arr) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(in) :: dim1 + real(c_double), intent(out) :: arr(dim1) + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + arr = c%lapse_rate +end subroutine + +subroutine adiabatclimate_lapse_rate_intended_get_size(ptr, dim1) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(out) :: dim1 + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + dim1 = size(c%lapse_rate_intended,1) +end subroutine + +subroutine adiabatclimate_lapse_rate_intended_get(ptr, dim1, arr) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(in) :: dim1 + real(c_double), intent(out) :: arr(dim1) + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + arr = c%lapse_rate_intended +end subroutine + +subroutine adiabatclimate_convective_newton_step_size_get(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + real(c_double), intent(out) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + val = c%convective_newton_step_size +end subroutine + +subroutine adiabatclimate_convective_newton_step_size_set(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + real(c_double), intent(in) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + c%convective_newton_step_size = val +end subroutine + subroutine adiabatclimate_rtol_get(ptr, val) bind(c) use clima, only: AdiabatClimate type(c_ptr), intent(in) :: ptr @@ -750,6 +869,96 @@ subroutine adiabatclimate_tol_make_column_set(ptr, val) bind(c) c%tol_make_column = val end subroutine +subroutine adiabatclimate_epsj_get(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + real(c_double), intent(out) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + val = c%epsj +end subroutine + +subroutine adiabatclimate_epsj_set(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + real(c_double), intent(in) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + c%epsj = val +end subroutine + +subroutine adiabatclimate_xtol_rc_get(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + real(c_double), intent(out) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + val = c%xtol_rc +end subroutine + +subroutine adiabatclimate_xtol_rc_set(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + real(c_double), intent(in) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + c%xtol_rc = val +end subroutine + +subroutine adiabatclimate_max_rc_iters_get(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(out) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + val = c%max_rc_iters +end subroutine + +subroutine adiabatclimate_max_rc_iters_set(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(in) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + c%max_rc_iters = val +end subroutine + +subroutine adiabatclimate_max_rc_iters_convection_get(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(out) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + val = c%max_rc_iters_convection +end subroutine + +subroutine adiabatclimate_max_rc_iters_convection_set(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + integer(c_int), intent(in) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + c%max_rc_iters_convection = val +end subroutine + +subroutine adiabatclimate_verbose_get(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + logical(c_bool), intent(out) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + val = c%verbose +end subroutine + +subroutine adiabatclimate_verbose_set(ptr, val) bind(c) + use clima, only: AdiabatClimate + type(c_ptr), intent(in) :: ptr + logical(c_bool), intent(in) :: val + type(AdiabatClimate), pointer :: c + call c_f_pointer(ptr, c) + c%verbose = val +end subroutine + subroutine adiabatclimate_p_surf_get(ptr, val) bind(c) use clima, only: AdiabatClimate type(c_ptr), intent(in) :: ptr diff --git a/src/adiabat/clima_adiabat.f90 b/src/adiabat/clima_adiabat.f90 index af5e3c8..b608de9 100644 --- a/src/adiabat/clima_adiabat.f90 +++ b/src/adiabat/clima_adiabat.f90 @@ -79,8 +79,8 @@ module clima_adiabat !> then the layer below is convecting with the current layer. Index 1 determines !> if the first atomspheric layer is convecting with the ground. logical, allocatable :: convecting_with_below(:) - real(dp), allocatable :: lapse_rate(:) !! The true lapse rate - real(dp), allocatable :: lapse_rate_intended(:) !! The computed lapse rate + real(dp), allocatable :: lapse_rate(:) !! The true lapse rate (dlnT/dlnP) + real(dp), allocatable :: lapse_rate_intended(:) !! The computed lapse rate (dlnT/dlnP) !> The size of the newton step. real(dp) :: convective_newton_step_size = 1.0e-1_dp @@ -95,8 +95,12 @@ module clima_adiabat real(dp) :: epsj = 1.0e-2_dp !> xtol for RC equilibrium real(dp) :: xtol_rc = 1.0e-5_dp + !> Max number of iterations in the RCE routine integer :: max_rc_iters = 10 - logical :: verbose = .true. + !> Max number of iterations for which convective layers can + !> be converged to radiative layers in the RCE routine + integer :: max_rc_iters_convection = 5 + logical :: verbose = .true. !! verbosity ! State of the atmosphere real(dp) :: P_surf !! Surface pressure (dynes/cm^2) @@ -152,14 +156,21 @@ module subroutine AdiabatClimate_make_profile_rc(self, P_i_surf, T_in, err) real(dp), intent(in) :: T_in(:) character(:), allocatable, intent(out) :: err end subroutine + + !> Compute full radiative-convective equilibrium. module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convecting_with_below, err) result(converged) class(AdiabatClimate), intent(inout) :: self + !> Array of surface pressures of each species (dynes/cm^2) real(dp), intent(in) :: P_i_surf(:) + !> A guess for the surface temperature (K) real(dp), intent(in) :: T_surf_guess + !> A guess for the temperature in each atmospheric layer (K) real(dp), intent(in) :: T_guess(:) + !> An array describing a guess for the radiative vs. convective + !> regions of the atmosphere logical, optional, intent(in) :: convecting_with_below(:) character(:), allocatable, intent(out) :: err - logical :: converged + logical :: converged !! Whether the routine converged or not. end function end interface diff --git a/src/adiabat/clima_adiabat_rc.f90 b/src/adiabat/clima_adiabat_rc.f90 index 4b1377b..36bc399 100644 --- a/src/adiabat/clima_adiabat_rc.f90 +++ b/src/adiabat/clima_adiabat_rc.f90 @@ -548,6 +548,16 @@ subroutine mixing_ratios(d, P, T, f_i_layer, f_dry) endif enddo + ! ng == 1 case + if (d%sp%ng == 1) then + f_i_layer(1) = 1.0_dp + if (d%sp_type(i) == CondensingSpeciesType) then + f_dry = 1.0e-40_dp + else + f_dry = 1.0_dp + endif + endif + end subroutine function dT_dP(d, P) result(dTdP) diff --git a/src/adiabat/clima_adiabat_solve.f90 b/src/adiabat/clima_adiabat_solve.f90 index 1e8286c..a8b33b7 100644 --- a/src/adiabat/clima_adiabat_solve.f90 +++ b/src/adiabat/clima_adiabat_solve.f90 @@ -85,9 +85,10 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec logical :: converged type(MinpackHybrj) :: mv - logical, allocatable :: convecting_with_below_save(:) + logical, allocatable :: convecting_with_below_save(:,:) + real(dp), allocatable :: difference(:) real(dp), allocatable :: T_in(:) - integer :: i + integer :: i, j if (.not.self%double_radiative_grid) then err = 'AdiabatClimate must be initialized with "double_radiative_grid" set to True '// & @@ -98,8 +99,9 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec err = "T_guess has the wrong dimension" return endif - - allocate(convecting_with_below_save(self%nz)) + + allocate(convecting_with_below_save(self%nz,0)) + allocate(difference(self%nz)) allocate(T_in(self%nz+1)) T_in(1) = T_surf_guess @@ -110,7 +112,7 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec call AdiabatClimate_set_convecting_zones(self, convecting_with_below, err) if (allocated(err)) return else - call AdiabatClimate_update_convecting_zones(self, P_i_surf, T_in, err) + call AdiabatClimate_update_convecting_zones(self, P_i_surf, T_in, .false., err) if (allocated(err)) return endif @@ -137,20 +139,54 @@ module function AdiabatClimate_RCE(self, P_i_surf, T_surf_guess, T_guess, convec err = 'hybrj root solve failed in surface_temperature: '//err return endif - convecting_with_below_save = self%convecting_with_below - call AdiabatClimate_update_convecting_zones(self, P_i_surf, mv%x, err) + + ! Update all variables to the current root + call AdiabatClimate_objective(self, P_i_surf, mv%x, .false., mv%fvec, err) if (allocated(err)) return - if (all(convecting_with_below_save .eqv. self%convecting_with_below)) then - ! No more convection changes, so converged + + ! Compute the difference between true and intended lapse rate + difference = self%lapse_rate - self%lapse_rate_intended + + ! Save the current convective zones + convecting_with_below_save = reshape(convecting_with_below_save,shape=[self%nz,i],pad=self%convecting_with_below) + + ! Update the convective zones + if (i < self%max_rc_iters_convection) then + ! Here, we permit + ! - convective layers to become radiative + ! - radiative layers to become convective + call AdiabatClimate_update_convecting_zones(self, P_i_surf, mv%x, .false., err) + if (allocated(err)) return + else + ! Here, we only permit radiative layers to become convective + call AdiabatClimate_update_convecting_zones(self, P_i_surf, mv%x, .true., err) + if (allocated(err)) return + endif + + ! Check for convergence + if (all(convecting_with_below_save(:,i) .eqv. self%convecting_with_below)) then + ! Convective zones did not change between iterations, so converged + converged = .true. + endif + if (any(difference > 1.0e-5_dp)) then + ! If there remains a layer that is superadiabatic (to within a tolerance) + ! then convergence has not been reached + converged = .false. + endif + if (converged) then if (self%verbose) then print'(1x,A)','CONVERGED' endif - converged = .true. exit endif enddo - + + ! Return convection information to what is was prior + ! to checking for convergence + call AdiabatClimate_set_convecting_zones(self, convecting_with_below_save(:,i-1), err) + if (allocated(err)) return + contains subroutine fcn(n_, x_, fvec_, fjac_, ldfjac_, iflag_) @@ -262,10 +298,10 @@ subroutine AdiabatClimate_objective_(self, P_i_surf, T_in, ignore_convection, re self%rad%f_total(1) = self%rad%f_total(1) + self%surface_heat_flow do i = 1,self%nz+1 - f_total(i) = self%rad%f_total(2*i-1) + f_total(i) = self%rad%f_total(2*i-1)/1.0e3_dp ! we normalized here for the purposes of optimization enddo - ! Radiative energy going into each layer (ergs/(cm^2*s)) + ! Radiative energy going into each layer fluxes(1) = f_total(1) do i = 2,self%nz+1 fluxes(i) = (f_total(i) - f_total(i-1)) @@ -276,16 +312,6 @@ subroutine AdiabatClimate_objective_(self, P_i_surf, T_in, ignore_convection, re return endif - ! j = self%ind_conv_upper(1) - ! ! Radiative layers - ! res(j+1:) = fluxes(j+1:) - ! ! Convective layers - ! i = j - 1 ! i is the atmospheric layer - ! res(j) = f_total(i+1) - ! do i = 1,j-1 - ! res(i) = self%lapse_rate(i) - self%lapse_rate_intended(i) - ! enddo - call AdiabatClimate_residuals_with_convection(self, f_total, self%lapse_rate, self%lapse_rate_intended, res) end subroutine @@ -392,11 +418,12 @@ subroutine AdiabatClimate_set_convecting_zones(self, convecting_with_below, err) !> tends to a T profile unstable to convection, then the layer is convective. The !> routine specifically updates self%convecting_with_below, self%n_convecting_zones, !> self%ind_conv_lower and self%ind_conv_upper. It also updates all atmosphere varibles. - subroutine AdiabatClimate_update_convecting_zones(self, P_i_surf, T_in, err) + subroutine AdiabatClimate_update_convecting_zones(self, P_i_surf, T_in, no_convection_to_radiation, err) use clima_useful, only: linear_solve class(AdiabatClimate), intent(inout) :: self real(dp), intent(in) :: P_i_surf(:) real(dp), intent(in) :: T_in(:) + logical, intent(in) :: no_convection_to_radiation character(:), allocatable, intent(out) :: err real(dp), allocatable :: F(:), dFdT(:,:), deltaT(:), T_perturb(:) @@ -433,13 +460,24 @@ subroutine AdiabatClimate_update_convecting_zones(self, P_i_surf, T_in, err) difference = lapse_rate_perturb - self%lapse_rate_intended - do i = 1,self%nz - if (difference(i) >= 0.0_dp) then - self%convecting_with_below(i) = .true. - else - self%convecting_with_below(i) = .false. - endif - enddo + if (no_convection_to_radiation) then + do i = 1,self%nz + if (.not.self%convecting_with_below(i)) then + difference(i) = self%lapse_rate(i) - self%lapse_rate_intended(i) + endif + if (difference(i) >= 0.0_dp) then + self%convecting_with_below(i) = .true. + endif + enddo + else + do i = 1,self%nz + if (difference(i) >= 0.0_dp) then + self%convecting_with_below(i) = .true. + else + self%convecting_with_below(i) = .false. + endif + enddo + endif call AdiabatClimate_set_convecting_zones(self, self%convecting_with_below, err) if (allocated(err)) return