surface albedo a vector. surface emissivity also included

clima/cython/Radtran.pyx

@@ -42,13 +42,34 @@ cdef class Radtran:
     return out_c[:-1].tobytes().decode()
 
   property surface_albedo:
-    "The surface albedo"
+    "ndarray[double,ndim=1]. The surface albedo in each solar wavelength bin"
     def __get__(self):
-      cdef double val
-      rad_pxd.radtran_surface_albedo_get(&self._ptr, &val)
-      return val
-    def __set__(self, double val):
-      rad_pxd.radtran_surface_albedo_set(&self._ptr, &val)
+      cdef int dim1
+      rad_pxd.radtran_surface_albedo_get_size(&self._ptr, &dim1)
+      cdef ndarray arr = np.empty(dim1, np.double)
+      rad_pxd.radtran_surface_albedo_get(&self._ptr, &dim1, <double *>arr.data)
+      return arr
+    def __set__(self, ndarray[double, ndim=1] arr):
+      cdef int dim1
+      rad_pxd.radtran_surface_albedo_get_size(&self._ptr, &dim1)
+      if arr.shape[0] != dim1:
+        raise ClimaException('"surface_albedo" is the wrong size')
+      rad_pxd.radtran_surface_albedo_set(&self._ptr, &dim1, <double *>arr.data)
+
+  property surface_emissivity:
+    "ndarray[double,ndim=1]. The surface emissivity in each IR wavelength bin"
+    def __get__(self):
+      cdef int dim1
+      rad_pxd.radtran_surface_emissivity_get_size(&self._ptr, &dim1)
+      cdef ndarray arr = np.empty(dim1, np.double)
+      rad_pxd.radtran_surface_emissivity_get(&self._ptr, &dim1, <double *>arr.data)
+      return arr
+    def __set__(self, ndarray[double, ndim=1] arr):
+      cdef int dim1
+      rad_pxd.radtran_surface_emissivity_get_size(&self._ptr, &dim1)
+      if arr.shape[0] != dim1:
+        raise ClimaException('"surface_emissivity" is the wrong size')
+      rad_pxd.radtran_surface_emissivity_set(&self._ptr, &dim1, <double *>arr.data)
 
   property ir:
     "The OpticalProperties for longwave radiative transfer"

clima/cython/Radtran_pxd.pxd

@@ -5,8 +5,13 @@ cdef extern void radtran_opacities2yaml_wrapper_1(void *ptr, int *out_len, void
 cdef extern void radtran_opacities2yaml_wrapper_2(void *ptr, void *out_cp, int *out_len, char* out_c)
 
 # getters and setters
-cdef extern void radtran_surface_albedo_get(void *ptr, double *val)
-cdef extern void radtran_surface_albedo_set(void *ptr, double *val)
+cdef extern void radtran_surface_albedo_get_size(void *ptr, int *dim1)
+cdef extern void radtran_surface_albedo_get(void *ptr, int *dim1, double *arr)
+cdef extern void radtran_surface_albedo_set(void *ptr, int *dim1, double *arr)
+
+cdef extern void radtran_surface_emissivity_get_size(void *ptr, int *dim1)
+cdef extern void radtran_surface_emissivity_get(void *ptr, int *dim1, double *arr)
+cdef extern void radtran_surface_emissivity_set(void *ptr, int *dim1, double *arr)
 
 cdef extern void radtran_ir_get(void *ptr, void *ptr1);
 

clima/fortran/Radtran.f90

@@ -50,22 +50,62 @@ subroutine radtran_opacities2yaml_wrapper_2(ptr, out_cp, out_len, out_c) bind(c)
 !!! getters and setters !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-subroutine radtran_surface_albedo_get(ptr, val) bind(c)
+subroutine radtran_surface_albedo_get_size(ptr, dim1) bind(c)
   use clima, only: Radtran
   type(c_ptr), intent(in) :: ptr
-  real(c_double), intent(out) :: val
+  integer(c_int), intent(out) :: dim1
   type(Radtran), pointer :: rad
   call c_f_pointer(ptr, rad)
-  val = rad%surface_albedo
+  dim1 = size(rad%surface_albedo,1)
 end subroutine
 
-subroutine radtran_surface_albedo_set(ptr, val) bind(c)
+subroutine radtran_surface_albedo_get(ptr, dim1, arr) bind(c)
   use clima, only: Radtran
   type(c_ptr), intent(in) :: ptr
-  real(c_double), intent(in) :: val
+  integer(c_int), intent(in) :: dim1
+  real(c_double), intent(out) :: arr(dim1)
   type(Radtran), pointer :: rad
   call c_f_pointer(ptr, rad)
-  rad%surface_albedo = val
+  arr = rad%surface_albedo
+end subroutine
+
+subroutine radtran_surface_albedo_set(ptr, dim1, arr) bind(c)
+  use clima, only: Radtran
+  type(c_ptr), intent(in) :: ptr
+  integer(c_int), intent(in) :: dim1
+  real(c_double), intent(in) :: arr(dim1)
+  type(Radtran), pointer :: rad
+  call c_f_pointer(ptr, rad)
+  rad%surface_albedo = arr
+end subroutine
+
+subroutine radtran_surface_emissivity_get_size(ptr, dim1) bind(c)
+  use clima, only: Radtran
+  type(c_ptr), intent(in) :: ptr
+  integer(c_int), intent(out) :: dim1
+  type(Radtran), pointer :: rad
+  call c_f_pointer(ptr, rad)
+  dim1 = size(rad%surface_emissivity,1)
+end subroutine
+
+subroutine radtran_surface_emissivity_get(ptr, dim1, arr) bind(c)
+  use clima, only: Radtran
+  type(c_ptr), intent(in) :: ptr
+  integer(c_int), intent(in) :: dim1
+  real(c_double), intent(out) :: arr(dim1)
+  type(Radtran), pointer :: rad
+  call c_f_pointer(ptr, rad)
+  arr = rad%surface_emissivity
+end subroutine
+
+subroutine radtran_surface_emissivity_set(ptr, dim1, arr) bind(c)
+  use clima, only: Radtran
+  type(c_ptr), intent(in) :: ptr
+  integer(c_int), intent(in) :: dim1
+  real(c_double), intent(in) :: arr(dim1)
+  type(Radtran), pointer :: rad
+  call c_f_pointer(ptr, rad)
+  rad%surface_emissivity = arr
 end subroutine
 
 subroutine radtran_ir_get(ptr, ptr1) bind(c)

examples/Tutorial_AdiabatClimate.ipynb

@@ -15,29 +15,18 @@
    "execution_count": 1,
    "id": "bfb3a831",
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "<threadpoolctl.threadpool_limits at 0x11e997190>"
-      ]
-     },
-     "execution_count": 1,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
     "import numpy as np\n",
     "import numba as nb\n",
     "from matplotlib import pyplot as plt\n",
     "from scipy import constants as const\n",
     "\n",
-    "# from clima import AdiabatClimate # If you installed clima alone from source\n",
-    "from photochem.clima import AdiabatClimate # If you installed photochem via conda\n",
+    "from clima import AdiabatClimate # If you installed clima alone from source\n",
+    "# from photochem.clima import AdiabatClimate # If you installed photochem via conda\n",
     "\n",
     "from threadpoolctl import threadpool_limits\n",
-    "threadpool_limits(limits=4) # set number of threads. Must come after import of AdiabatClimate"
+    "_ = threadpool_limits(limits=4) # set number of threads. Must come after import of AdiabatClimate"
    ]
   },
   {
@@ -63,7 +52,8 @@
     "\n",
     "c = AdiabatClimate(species_file, \n",
     "                   settings_file, \n",
-    "                   star_file)"
+    "                   star_file,\n",
+    "                   double_radiative_grid = False)"
    ]
   },
   {
@@ -629,7 +619,7 @@
    "outputs": [],
    "source": [
     "c.T_trop = 157.0 # Follows Kopparapu et al.\n",
-    "c.rad.surface_albedo = 0.2\n",
+    "c.rad.surface_albedo = np.ones(c.rad.surface_albedo.shape[0])*0.2\n",
     "\n",
     "P_CO2s = np.logspace(0,2,100)\n",
     "P_i = np.array([200, 8, 1.0, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10])*1e6\n",
@@ -725,7 +715,7 @@
    "outputs": [],
    "source": [
     "# return values to what they were before\n",
-    "c.rad.surface_albedo = 0.25\n",
+    "c.rad.surface_albedo = np.ones(c.rad.surface_albedo.shape[0])*0.25\n",
     "c.T_trop = 215.0"
    ]
   },
@@ -863,7 +853,7 @@
     "    ISR_, OLR_ = c.TOA_fluxes_column(T, N_i)\n",
     "    OLR[i] = OLR_/1.0e3\n",
     "    ISR[i] = ISR_/1.0e3\n",
-    "    albedo[i] = c.rad.surface_albedo"
+    "    albedo[i] = c.rad.surface_albedo[0]"
    ]
   },
   {
@@ -924,7 +914,7 @@
    "source": [
     "c.T_trop = 215 # restore T_trop\n",
     "c.albedo_fcn = None # remove albedo feedback\n",
-    "c.rad.surface_albedo = 0.25 # restore albedo"
+    "c.rad.surface_albedo = np.ones(c.rad.surface_albedo.shape[0])*0.25 # restore albedo"
    ]
   },
   {

src/clima.f90

@@ -10,7 +10,7 @@ module clima
   use clima_eqns, only: temp_dependent_albedo_fcn, ocean_solubility_fcn
 
   ! Radiative transfer classes
-  use clima_radtran, only: Radtran, RadtranIR
+  use clima_radtran, only: Radtran
 
   ! Adiabatic climate model
   use clima_adiabat, only: AdiabatClimate