rnatm

commit 3053c964f897602e875acbe6a9f2d173f3378145
parent 24325de981241a53d0d3748f6141a0a9c1202cbb
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 30 Aug 2022 10:41:57 +0200

Fix the calculation of the radiative coefficients of aerosols

Diffstat:
M
src/rnatm_octree.c
 | 
58
++++++++++++++++++++++++++++++++++++++++------------------

1 file changed, 40 insertions(+), 18 deletions(-)
diff --git a/src/rnatm_octree.c b/src/rnatm_octree.c
@@ -500,7 +500,6 @@ setup_tetra_radcoefs_aerosol
   (struct rnatm* atm,
    const struct suvm_primitive* tetra,
    const size_t iband,
-   const size_t iquad_pt,
    const struct aerosol* aerosol,
    struct radcoefs* radcoefs)
 {
@@ -508,6 +507,7 @@ setup_tetra_radcoefs_aerosol
   struct sars_band ars_band;
   double gas_spectral_range[2]; /* In cm^-1 */
   size_t ars_ibands[2];
+  float ka[4], ks[4], kext[4];
   ASSERT(atm && aerosol && tetra && radcoefs);
 
   /* Look for the aerosol bands covered by the gas band */
@@ -516,35 +516,56 @@ setup_tetra_radcoefs_aerosol
   gas_spectral_range[1] = nextafter(gas_band.upper, 0); /* inclusive */
   SARS(find_bands(aerosol->sars, gas_spectral_range, ars_ibands));
 
-  /* TODO ask VE what to do if:
-   *  - the spectral meshes have holes
-   *  - no aerosol band is overlaid by the gas band => FIXME no aerosol
-   *  - the aerosol bands are entirely included in the gas band FIXME general case
-   *  - gas degenerated => FIXME the aerosol bands must include the gas band
-   *  - the gas band and aerosol band are the same => FIXME take the aerosol properties */
-
-  /* Aerosol spectral data is not overlaid by the gas band */
-  if(ars_ibands[0] > ars_ibands[1]) FATAL("Not implemented yet!\n");
+  /* No aerosol band is overlaid by the gas band: ignore the aerosol */
+  if(ars_ibands[0] > ars_ibands[1]) {
+    radcoefs->ks_min = FLT_MAX;
+    radcoefs->ks_max =-FLT_MAX;
+    radcoefs->ka_min = FLT_MAX;
+    radcoefs->ka_max =-FLT_MAX;
+    radcoefs->kext_min = FLT_MAX;
+    radcoefs->kext_max =-FLT_MAX;
+    return;
+  }
 
   SARS(get_band(aerosol->sars, ars_ibands[0], &ars_band));
 
-  /* The gas band is included in an aerosol band: use gas spectral data */
+  /* The gas band is included in an aerosol band: use the aerosol radiative
+   * coefficients */
   if(ars_ibands[0] == ars_ibands[1]
   && ars_band.lower <= gas_band.lower
   && ars_band.upper >= gas_band.upper) {
-    /* FIXME use the aerosole properties */
-    setup_tetra_radcoefs_gas(atm, tetra, iband, iquad_pt, radcoefs);
+
+    /* Compute the scattering coefficient range of the tetrahedron */
+    ks[0] = ars_band.ks_list[tetra->indices[0]];
+    ks[1] = ars_band.ks_list[tetra->indices[1]];
+    ks[2] = ars_band.ks_list[tetra->indices[2]];
+    ks[3] = ars_band.ks_list[tetra->indices[3]];
+    radcoefs->ks_min = MMIN(MMIN(ks[0], ks[1]), MMIN(ks[2], ks[3]));
+    radcoefs->ks_max = MMAX(MMAX(ks[0], ks[1]), MMAX(ks[2], ks[3]));
+
+    /* Compute the absorption coefficient range of the tetrahedron */
+    ka[0] = ars_band.ka_list[tetra->indices[0]];
+    ka[1] = ars_band.ka_list[tetra->indices[1]];
+    ka[2] = ars_band.ka_list[tetra->indices[2]];
+    ka[3] = ars_band.ka_list[tetra->indices[3]];
+    radcoefs->ka_min = MMIN(MMIN(ka[0], ka[1]), MMIN(ka[2], ka[3]));
+    radcoefs->ka_max = MMAX(MMAX(ka[0], ka[1]), MMAX(ka[2], ka[3]));
+
+    /* Compute the extinction coefficient range of the tetrahedron */
+    kext[0] = ka[0] + ks[0];
+    kext[1] = ka[1] + ks[1];
+    kext[2] = ka[2] + ks[2];
+    kext[3] = ka[3] + ks[3];
+    radcoefs->kext_min = MMIN(MMIN(kext[0], kext[1]), MMIN(kext[2], kext[3]));
+    radcoefs->kext_max = MMAX(MMAX(kext[0], kext[1]), MMAX(kext[2], kext[3]));
 
   /* The gas band overlaid N aerosol bands (N >= 1) */
   } else {
     float tau_ka[4] = {0, 0, 0, 0};
     float tau_ks[4] = {0, 0, 0, 0};
-    float ka[4], ks[4], kext[4];
     float rcp_gas_band_len;
     size_t iars_band;
 
-    /* TODO If gas.lambda_min == gas.lambda_max use the properties of the
-     * aerosol band */
     FOR_EACH(iars_band, ars_ibands[0], ars_ibands[1]+1) {
       double lambda_min;
       double lambda_max;
@@ -554,6 +575,7 @@ setup_tetra_radcoefs_aerosol
       lambda_min = MMAX(gas_band.lower, ars_band.lower);
       lambda_max = MMIN(gas_band.upper, ars_band.upper);
       lambda_len = (float)(lambda_max - lambda_min);
+      ASSERT(lambda_len > 0);
 
       tau_ks[0] += ars_band.ks_list[tetra->indices[0]] * lambda_len;
       tau_ks[1] += ars_band.ks_list[tetra->indices[1]] * lambda_len;
@@ -567,7 +589,7 @@ setup_tetra_radcoefs_aerosol
     }
 
     /* Compute the radiative coefficients of the tetrahedron */
-    ASSERT(gas_band.upper != gas_band.lower); /* FIXME is it possible? */
+    ASSERT(gas_band.upper > gas_band.lower);
     rcp_gas_band_len = 1.f / (float)(gas_band.upper - gas_band.lower);
     f4_mulf(ks, tau_ks, rcp_gas_band_len);
     f4_mulf(ka, tau_ka, rcp_gas_band_len);
@@ -599,7 +621,7 @@ setup_tetra_radcoefs
     setup_tetra_radcoefs_gas(atm, tetra, iband, iquad_pt, radcoefs);
   } else {
     const struct aerosol* aerosol = darray_aerosol_cdata_get(&atm->aerosols)+cpnt;
-    setup_tetra_radcoefs_aerosol(atm, tetra, iband, iquad_pt, aerosol, radcoefs);
+    setup_tetra_radcoefs_aerosol(atm, tetra, iband, aerosol, radcoefs);
   }
 }