rnatm

commit c32e740b4dcbcea535824e2f41bc3a464b11efde
parent fb6684e99fccc01a7487ed78ac845d11ca71283c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sat,  1 Oct 2022 17:56:55 +0200

Rm aerosols with no radiative props for input spectral domain

Diffstat:
M
src/rnatm_octrees_storage.c
 | 
11
++++++-----
M
src/rnatm_properties.c
 | 
92
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-----------------

2 files changed, 79 insertions(+), 24 deletions(-)
diff --git a/src/rnatm_octrees_storage.c b/src/rnatm_octrees_storage.c
@@ -152,11 +152,12 @@ compute_octrees_signature(const struct rnatm* atm, hash256_T signature)
       res = sars_find_bands(aerosol->sars, atm->spectral_range, ibands);
       if(res != RES_OK) goto error;
 
-      if(ibands[0] <= ibands[1]) {
-        /* The aerosol is overlapped by the spectral range to consider */
-        res = compute_aerosol_signature(aerosol, ibands, aerosol_signatures[i]);
-        if(res != RES_OK) goto error;
-      }
+      /* Note that if an aerosol had no data for the input spectral range, it
+       * was suppressed when setting the atmosphere properties */
+      ASSERT(ibands[0] <= ibands[1]);
+
+      res = compute_aerosol_signature(aerosol, ibands, aerosol_signatures[i]);
+      if(res != RES_OK) goto error;
     }
 
     /* Pay attention to the order of how aerosol hashes are registered.
diff --git a/src/rnatm_properties.c b/src/rnatm_properties.c
@@ -588,7 +588,7 @@ setup_gas_properties(struct rnatm* atm, const struct rnatm_gas_args* gas_args)
   if(res != RES_OK) goto error;
   res = sck_get_band(atm->gas.ck, nbands-1, &band_upp);
   if(res != RES_OK) goto error;
-  log_info(atm, "the gas is composed of %lu band%sin [%g, %g] nm\n",
+  log_info(atm, "the gas is composed of %lu band%sin [%g, %g[ nm\n",
     nbands, nbands > 1 ? "s " : " ",
     band_low.lower,
     band_upp.upper);
@@ -607,6 +607,7 @@ setup_aerosol_properties
    const struct rnatm_aerosol_args* aerosol_args)
 {
   char buf[128];
+  size_t ibands[2];
   struct time t0, t1;
   struct sars_create_args sars_args = SARS_CREATE_ARGS_DEFAULT;
   struct sbuf_create_args sbuf_args = SBUF_CREATE_ARGS_DEFAULT;
@@ -619,6 +620,28 @@ setup_aerosol_properties
   log_info(atm, "load %s properties\n", str_cget(&aerosol->name));
   time_current(&t0);
 
+  /* Create the Star-Aerosol loader */
+  sars_args.logger = atm->logger;
+  sars_args.allocator = atm->allocator;
+  sars_args.verbose = atm->verbose;
+  res = sars_create(&sars_args, &aerosol->sars);
+  if(res != RES_OK) goto error;
+
+  /* Load the aerosol radiative properties */
+  res = sars_load(aerosol->sars, aerosol_args->sars_filename);
+  if(res != RES_OK) goto error;
+  res = check_aerosol_sars_desc(atm, aerosol, aerosol_args);
+  if(res != RES_OK) goto error;
+
+  /* Ignore the aerosol if it has no data for the input spectral range */
+  res = sars_find_bands(aerosol->sars, atm->spectral_range, ibands);
+  if(res != RES_OK) goto error;
+  if(ibands[0] > ibands[1]) {
+    reset_aerosol_properties(aerosol);
+    goto exit;
+  }
+
+  /* Load the aerosol phase functions */
   res = load_phase_fn_list(atm, aerosol, aerosol_args);
   if(res != RES_OK) goto error;
 
@@ -637,23 +660,6 @@ setup_aerosol_properties
   res = check_aerosol_phase_fn_ids_desc(atm, aerosol, &sbuf_desc, aerosol_args);
   if(res != RES_OK) goto error;
 
-  /* Create the Star-Aerosol loader */
-  sars_args.logger = atm->logger;
-  sars_args.allocator = atm->allocator;
-  sars_args.verbose = atm->verbose;
-  res = sars_create(&sars_args, &aerosol->sars);
-  if(res != RES_OK) goto error;
-
-  /* Load the aerosol radiative properties */
-  res = sars_load(aerosol->sars, aerosol_args->sars_filename);
-  if(res != RES_OK) goto error;
-  res = check_aerosol_sars_desc(atm, aerosol, aerosol_args);
-  if(res != RES_OK) goto error;
-
-  /* TODO remove the aerosol if its spectral data are not covered by the
-   * spectral domain, i.e. the aerosol has no influence on the radiative
-   * transfer for the input spectral domain */
-
   /* Print elapsed time */
   time_sub(&t0, time_current(&t1), &t0);
   time_dump(&t0, TIME_ALL, NULL, buf, sizeof(buf));
@@ -668,6 +674,50 @@ error:
 }
 
 static res_T
+remove_useless_aerosols(struct rnatm* atm)
+{
+  struct aerosol* aerosols = NULL;
+  size_t i, n;
+  res_T res = RES_OK;
+  ASSERT(atm);
+
+  aerosols = darray_aerosol_data_get(&atm->aerosols);
+
+  n = 0;
+  FOR_EACH(i, 0, darray_aerosol_size_get(&atm->aerosols)) {
+    /* Discard aerosols with no radiative properties */
+    if(!aerosols[i].sars) {
+      log_warn(atm,
+        "discard %s aerosol: no radiative properties are defined for the "
+        "spectral domain [%g, %g] nm\n",
+        str_cget(&aerosols[i].name),
+        atm->spectral_range[0],
+        atm->spectral_range[1]);
+      continue;
+    }
+    /* Compact the aerosols to consider */
+    if(i != n) {
+      aerosol_copy_and_release(aerosols+n, aerosols+i);
+      ++n;
+    }
+  }
+
+  if(!n) {
+    /* Discard all aerosols */
+    darray_aerosol_purge(&atm->aerosols);
+  } else {
+    /* Resize the aerosols list */
+    res = darray_aerosol_resize(&atm->aerosols, n);
+    if(res != RES_OK) goto error;
+  }
+
+exit:
+  return res;
+error:
+  goto exit;
+}
+
+static res_T
 find_band_range
   (const struct rnatm* atm,
    const double range[2], /* In nanometers */
@@ -893,13 +943,17 @@ setup_properties(struct rnatm* atm, const struct rnatm_create_args* args)
   res = setup_gas_properties(atm, &args->gas);
   if(res != RES_OK) goto error;
 
+  res = setup_spectral_range(atm, args);
+  if(res != RES_OK) goto error;
+
   FOR_EACH(i, 0, darray_aerosol_size_get(&atm->aerosols)) {
     struct aerosol* aerosol = darray_aerosol_data_get(&atm->aerosols)+i;
     res = setup_aerosol_properties(atm, aerosol, args->aerosols+i);
     if(res != RES_OK) goto error;
   }
 
-  res = setup_spectral_range(atm, args);
+  /* Remove aerosols with no radiative properties for the input spectral range */
+  res = remove_useless_aerosols(atm);
   if(res != RES_OK) goto error;
 
 exit: