/*
   Name- convertFromLevelsToLayers.c

   Language- C     Type- MAIN

   Version- 1.0    Date-  2/28/2026   Programmer- Mike Pettey (STC)

   Function- This function interpolates a data value to a given pressure.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void interpolateProfile(int num_before_levels, float *before_pressures, float *before_data,
                        int num_after_levels, float *after_pressures, float *after_data);


void convertFromLevelsToLayers(int num_before_levels, float *before_pressures, float *before_data,
                               int num_after_layers, float *after_pressures, float *after_data)
  {
  int    layer, level, index, num_interpolated_levels, num_levels, nlayer;
  float  top_pressure, bottom_pressure, surface_pressure, surface_data, sumwgt, wgt;
  float  interpolated_data[100];
  float  *fine_layer_means;

  float interpolated_pressures[] = {  0.005, 0.025, 0.055, 0.103, 0.176, 0.280, 0.420, 0.604, 0.837,
                                      1.128, 1.484, 1.910, 2.416, 3.008, 3.696, 4.485, 5.384, 6.401,
                                      7.544, 8.821, 10.238, 11.782, 13.480, 15.378, 17.429, 19.676,
                                     22.127, 24.725, 27.573, 30.624, 33.872, 37.371, 41.119, 45.070,
                                     49.267, 53.767, 58.467, 63.465, 68.765, 74.314, 80.163, 86.312,
                                     92.761, 99.510, 106.559, 113.909, 121.607, 129.657, 138.006,
                                     146.705, 155.755, 165.153, 174.952, 185.102, 195.552, 206.400,
                                     217.650, 229.249, 241.249, 253.599, 266.347, 279.497, 293.046,
                                     306.996, 321.295, 336.044, 351.195, 366.743, 382.693, 399.042,
                                     415.842, 433.042, 450.691, 468.742, 487.191, 506.090, 525.389,
                                     545.089, 565.238, 585.838, 606.838, 628.238, 650.039, 672.288,
                                     694.988, 718.086, 741.637, 765.586, 789.936, 814.735, 839.935,
                                     865.586, 891.635, 918.085, 944.936, 972.184, 1000.000, 1027.987,
                                    1056.486, 1085.384 };

  num_interpolated_levels = 100;

  // Initialize the output array with missing values

  for (layer=0; layer<num_after_layers; layer++)
    after_data[layer] = 32768.0;


  // Determine the bottom pressure of the data which is defined as the lowest pressure
  // level containing non-missing data

  surface_pressure = 0.0;

  for (level=0; level<num_before_levels; level++)
    {
    if (before_data[level] != -32768.0)
      surface_pressure = before_pressures[level];
    }

  if (surface_pressure > 1100.0)
    surface_pressure = 1100.0;


  // Interpolate the incoming level data to the 100 levels

  interpolateProfile(num_before_levels, before_pressures, before_data,
                     num_interpolated_levels, interpolated_pressures, interpolated_data);


  // Calculate the fine layer mean data

  fine_layer_means = (float*)malloc(num_interpolated_levels*sizeof(float));

  for (level=0; level<num_interpolated_levels; level++)
    fine_layer_means[level] = -32768.0;

  num_levels = num_interpolated_levels;

  if (interpolated_data[num_levels-1] == -32768.0)
    {
    index = num_interpolated_levels - 1;

    while ((index >= 0) && (interpolated_data[index] == -32768.0))
      {
      index--;
      num_levels--;
      }
    }

  if (num_levels > 0)
    {
    surface_data = -32768.0;

    if ((interpolated_data[num_levels-2] != -32768.0) && (interpolated_data[num_levels-1] != -32768.0))
      {
      surface_data = interpolated_data[num_levels-2] +
                     (interpolated_data[num_levels-1] - interpolated_data[num_levels-2]) *
                     (surface_pressure - interpolated_pressures[num_levels-2]) /
                     (interpolated_pressures[num_levels-1]-before_pressures[num_levels-2]);
      }

    fine_layer_means[0] = interpolated_data[0];

    for (level=1; level<(num_levels-1); level++)
      {
      if ((interpolated_data[level] != -32768.0) && (interpolated_data[level-1] != -32768.0))
        fine_layer_means[level] = 0.5 * (interpolated_data[level-1] + interpolated_data[level]);
      }

    if ((interpolated_data[num_levels-1] != -32768.0) && (surface_data != -32768.0))
      fine_layer_means[num_levels-1] = 0.5 * (interpolated_data[num_levels-1] + surface_data);
    }


  // Calculate the slab mean (layer) values

  for (layer=0; layer<num_after_layers; layer++)
    {
    top_pressure    = after_pressures[layer];
    bottom_pressure = after_pressures[layer+1];

    after_data[layer] = 0.0;
    nlayer = 0;
    sumwgt = 0.0;

    // Weight by the thickness of layers to account for the fact that
    // the layer thickness varies in the pressure grid

    for (level=1; level<num_interpolated_levels; level++)
      {
      if (interpolated_pressures[level] <= surface_pressure)
        {
        if ((interpolated_pressures[level] > top_pressure) && (interpolated_pressures[level] <= bottom_pressure))
          {
          if (interpolated_data[level] != -32768.0)
            {
            nlayer++;
            wgt = log(interpolated_pressures[level] / interpolated_pressures[level-1]);
            after_data[layer] += (wgt * interpolated_data[level]);

            sumwgt += wgt;
            }
          }
        }
      }

    if (nlayer > 0)
      after_data[layer] = after_data[layer] / sumwgt;
    else
      after_data[layer] = -32768.0;
    }

  free(fine_layer_means);
  }

// end of file