added dry routines to make atmosphere

clima/cython/AdiabatClimate.pyx

@@ -87,6 +87,37 @@ cdef class AdiabatClimate:
     &ng, <double *>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, <double *>P.data,
+      &dim_T, <double *>T.data,
+      &dim1_f_i, &dim2_f_i, <double *>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, <double *>P.data,
+      &dim_T, <double *>T.data,
+      &dim1_f_i, &dim2_f_i, <double *>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

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)
 

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

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
 )
 

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`.