<commit_before>#include "snow_drift.h"
#include "forecast_time.h"
#include "level.h"
#include "logger.h"
#include "metutil.h"
#include "util.h"

#include "fetcher.h"
#include "hitool.h"
#include "radon.h"
#include "writer.h"

using namespace himan;
using namespace himan::plugin;

const param SDIParam("SNOWDRIFT-N");
const param SAParam("SNACC-H");
const param DAParam("SNDACC-N");
const double SFLimit = 0.09;

double DriftMagnitude(double ff, double mi);
double DriftIndex(double sv);
double MobilityIndex(double da, double sa);

void CalculateSnowDriftIndex(info_t& myTargetInfo, const std::vector<double>& T, const std::vector<double>& FF,
                             const std::vector<double>& SF, const std::vector<double>& pSA,
                             const std::vector<double>& pDA)
{
	myTargetInfo->Find<param>(SDIParam);
	auto& SI = VEC(myTargetInfo);
	myTargetInfo->Find<param>(SAParam);
	auto& SA = VEC(myTargetInfo);
	myTargetInfo->Find<param>(DAParam);
	auto& DA = VEC(myTargetInfo);

	for (auto&& tup : zip_range(SI, SA, DA, T, FF, SF, pSA, pDA))
	{
		auto& si = tup.get<0>();        // snow drift index
		auto& sa = tup.get<1>();        // accumulated snow cover age
		auto& da = tup.get<2>();        // accumulated snowdrift value
		const auto t = tup.get<3>();    // temperature
		const auto ff = tup.get<4>();   // wind speed
		const auto sf = tup.get<5>();   // snow fall rate (mm/h)
		const auto psa = tup.get<6>();  // previous accumulated snow cover age
		const auto pda = tup.get<7>();  // previous accumulated snowdrift value

		if (IsMissing(t) || IsMissing(ff) || IsMissing(sf))
		{
			continue;
		}

		if (t > 273.15)
		{
			// no possibility for snow drift until next snow fall
			si = 0;
			sa = 0;
			da = 0;
			continue;
		}

		/*
		        this is not present in the fortran code either
		        if (ff < 6.0)
		        {
		            si = 0;
		        }
		*/

		if (sf > SFLimit)
		{
			// during snow fall
			si = DriftIndex(DriftMagnitude(ff, 1.0));
			sa = 0;
			da = 0;
		}
		else
		{
			// after snow fall
			sa = psa + 1;  // add one hour to snowcover age

			if (ff < 6)
			{
				si = 0;
				da = pda;
			}
			else
			{
				const double sv = DriftMagnitude(ff, MobilityIndex(pda, sa));
				da = pda + sv;
				si = DriftIndex(sv);
			}
		}
	}
}

snow_drift::snow_drift()
{
	itsLogger = logger("snow_drift");
}

void snow_drift::Process(std::shared_ptr<const plugin_configuration> conf)
{
	Init(conf);

	SetParams({SDIParam, SAParam, DAParam});

	itsThreadDistribution = ThreadDistribution::kThreadForForecastTypeAndLevel;

	Start();
}

void snow_drift::Calculate(std::shared_ptr<info<double>> myTargetInfo, unsigned short theThreadIndex)
{
	auto myThreadedLogger = logger("snow_driftThread #" + std::to_string(theThreadIndex));

	const auto forecastTime = myTargetInfo->Time();
	const auto forecastLevel = myTargetInfo->Level();
	const auto forecastType = myTargetInfo->ForecastType();
	const auto prod = myTargetInfo->Producer();

	if (forecastTime.Step() == 0 && prod.Id() != 109)
	{
		// We only calculate analysis hour for LAPS
		return;
	}

	if (forecastTime.StepResolution() != kHourResolution || itsConfiguration->ForecastStep() != 1)
	{
		myThreadedLogger.Error("Snow drift can only be calculated with one hour resolution");
		return;
	}

	myThreadedLogger.Info("Calculating time " + static_cast<std::string>(forecastTime.ValidDateTime()) + " level " +
	                      static_cast<std::string>(forecastLevel));

	const std::string deviceType = "CPU";

	auto TInfo = Fetch(forecastTime, level(kHeight, 2), param("T-K"), forecastType, false);
	auto FFInfo = Fetch(forecastTime, level(kHeight, 10), param("FF-MS"), forecastType, false);
	auto SFInfo = Fetch(forecastTime, level(kHeight, 0), param("SNR-KGM2"), forecastType, false);

	if (!TInfo || !FFInfo || !SFInfo)
	{
		return;
	}

	auto prevTime = forecastTime;

	std::shared_ptr<info<double>> pSAInfo, pDAInfo;

	prevTime.ValidDateTime().Adjust(kHourResolution, -1);

	// In the start of the forecast fetch the latest sa and da
	// values from LAPS producer.

	if (prod.Id() == 109 || forecastTime.Step() == 1)
	{
		prevTime.OriginDateTime(prevTime.ValidDateTime());

		// Allow up to three hours of missing data
		for (int i = 0; i < 3; i++)
		{
			pSAInfo = Fetch(prevTime, level(kHeight, 0), SAParam, forecastType, itsConfiguration->SourceGeomNames(),
			                producer(109), false);
			pDAInfo = Fetch(prevTime, level(kHeight, 0), DAParam, forecastType, itsConfiguration->SourceGeomNames(),
			                producer(109), false);

			if (pSAInfo && pDAInfo)
			{
				break;
			}

			prevTime.ValidDateTime().Adjust(kHourResolution, -1);
			prevTime.OriginDateTime(prevTime.ValidDateTime());
		}
	}
	else
	{
		pSAInfo = Fetch(prevTime, level(kHeight, 0), SAParam, forecastType, false);
		pDAInfo = Fetch(prevTime, level(kHeight, 0), DAParam, forecastType, false);
	}

	if (!pSAInfo || !pDAInfo)
	{
		if (prod.Id() != 109)
		{
			myThreadedLogger.Error("DA and SA not found from obs producer");
			return;<commit_msg>Specify producer name with the id<commit_after>#include "snow_drift.h"
#include "forecast_time.h"
#include "level.h"
#include "logger.h"
#include "metutil.h"
#include "util.h"

#include "fetcher.h"
#include "hitool.h"
#include "radon.h"
#include "writer.h"

using namespace himan;
using namespace himan::plugin;

const param SDIParam("SNOWDRIFT-N");
const param SAParam("SNACC-H");
const param DAParam("SNDACC-N");
const double SFLimit = 0.09;

double DriftMagnitude(double ff, double mi);
double DriftIndex(double sv);
double MobilityIndex(double da, double sa);

void CalculateSnowDriftIndex(info_t& myTargetInfo, const std::vector<double>& T, const std::vector<double>& FF,
                             const std::vector<double>& SF, const std::vector<double>& pSA,
                             const std::vector<double>& pDA)
{
	myTargetInfo->Find<param>(SDIParam);
	auto& SI = VEC(myTargetInfo);
	myTargetInfo->Find<param>(SAParam);
	auto& SA = VEC(myTargetInfo);
	myTargetInfo->Find<param>(DAParam);
	auto& DA = VEC(myTargetInfo);

	for (auto&& tup : zip_range(SI, SA, DA, T, FF, SF, pSA, pDA))
	{
		auto& si = tup.get<0>();        // snow drift index
		auto& sa = tup.get<1>();        // accumulated snow cover age
		auto& da = tup.get<2>();        // accumulated snowdrift value
		const auto t = tup.get<3>();    // temperature
		const auto ff = tup.get<4>();   // wind speed
		const auto sf = tup.get<5>();   // snow fall rate (mm/h)
		const auto psa = tup.get<6>();  // previous accumulated snow cover age
		const auto pda = tup.get<7>();  // previous accumulated snowdrift value

		if (IsMissing(t) || IsMissing(ff) || IsMissing(sf))
		{
			continue;
		}

		if (t > 273.15)
		{
			// no possibility for snow drift until next snow fall
			si = 0;
			sa = 0;
			da = 0;
			continue;
		}

		/*
		        this is not present in the fortran code either
		        if (ff < 6.0)
		        {
		            si = 0;
		        }
		*/

		if (sf > SFLimit)
		{
			// during snow fall
			si = DriftIndex(DriftMagnitude(ff, 1.0));
			sa = 0;
			da = 0;
		}
		else
		{
			// after snow fall
			sa = psa + 1;  // add one hour to snowcover age

			if (ff < 6)
			{
				si = 0;
				da = pda;
			}
			else
			{
				const double sv = DriftMagnitude(ff, MobilityIndex(pda, sa));
				da = pda + sv;
				si = DriftIndex(sv);
			}
		}
	}
}

snow_drift::snow_drift()
{
	itsLogger = logger("snow_drift");
}

void snow_drift::Process(std::shared_ptr<const plugin_configuration> conf)
{
	Init(conf);

	SetParams({SDIParam, SAParam, DAParam});

	itsThreadDistribution = ThreadDistribution::kThreadForForecastTypeAndLevel;

	Start();
}

void snow_drift::Calculate(std::shared_ptr<info<double>> myTargetInfo, unsigned short theThreadIndex)
{
	auto myThreadedLogger = logger("snow_driftThread #" + std::to_string(theThreadIndex));

	const auto forecastTime = myTargetInfo->Time();
	const auto forecastLevel = myTargetInfo->Level();
	const auto forecastType = myTargetInfo->ForecastType();
	const auto prod = myTargetInfo->Producer();

	if (forecastTime.Step() == 0 && prod.Id() != 109)
	{
		// We only calculate analysis hour for LAPS
		return;
	}

	if (forecastTime.StepResolution() != kHourResolution || itsConfiguration->ForecastStep() != 1)
	{
		myThreadedLogger.Error("Snow drift can only be calculated with one hour resolution");
		return;
	}

	myThreadedLogger.Info("Calculating time " + static_cast<std::string>(forecastTime.ValidDateTime()) + " level " +
	                      static_cast<std::string>(forecastLevel));

	const std::string deviceType = "CPU";

	auto TInfo = Fetch(forecastTime, level(kHeight, 2), param("T-K"), forecastType, false);
	auto FFInfo = Fetch(forecastTime, level(kHeight, 10), param("FF-MS"), forecastType, false);
	auto SFInfo = Fetch(forecastTime, level(kHeight, 0), param("SNR-KGM2"), forecastType, false);

	if (!TInfo || !FFInfo || !SFInfo)
	{
		return;
	}

	auto prevTime = forecastTime;

	std::shared_ptr<info<double>> pSAInfo, pDAInfo;

	prevTime.ValidDateTime().Adjust(kHourResolution, -1);

	// In the start of the forecast fetch the latest sa and da
	// values from LAPS producer.

	if (prod.Id() == 109 || forecastTime.Step() == 1)
	{
		prevTime.OriginDateTime(prevTime.ValidDateTime());

		// Allow up to three hours of missing data
		for (int i = 0; i < 3; i++)
		{
			pSAInfo = Fetch(prevTime, level(kHeight, 0), SAParam, forecastType, itsConfiguration->SourceGeomNames(),
			                producer(109, 0, 0, "LAPSSCAN"), false);
			pDAInfo = Fetch(prevTime, level(kHeight, 0), DAParam, forecastType, itsConfiguration->SourceGeomNames(),
			                producer(109, 0, 0, "LAPSSCAN"), false);

			if (pSAInfo && pDAInfo)
			{
				break;
			}

			prevTime.ValidDateTime().Adjust(kHourResolution, -1);
			prevTime.OriginDateTime(prevTime.ValidDateTime());
		}
	}
	else
	{
		pSAInfo = Fetch(prevTime, level(kHeight, 0), SAParam, forecastType, false);
		pDAInfo = Fetch(prevTime, level(kHeight, 0), DAParam, forecastType, false);
	}

	if (!pSAInfo || !pDAInfo)
	{
		if (prod.Id() != 109)
		{
			myThreadedLogger.Error("DA and SA not found from obs producer");
			return;
{'fraction_non_alphanumeric': 0.14440504334121354, 'fraction_numerical': 0.012214342001576044, 'mean_word_length': 3.5004432624113475, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2213957', 'n_tokens_mistral': 3931, 'n_tokens_neox': 3483, 'n_words': 886}