Skip to content

Commit

Permalink
Code revisions.
Browse files Browse the repository at this point in the history
  • Loading branch information
@lars2015
lars2015 committed Jan 16, 2025
1 parent 6460415 commit a25a378
Show file tree
Hide file tree
Showing 30 changed files with 225 additions and 188 deletions.
2 changes: 1 addition & 1 deletion src/atm2grid.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
You should have received a copy of the GNU General Public License
along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2013-2024 Forschungszentrum Juelich GmbH
Copyright (C) 2013-2025 Forschungszentrum Juelich GmbH
*/

/*!
Expand Down
2 changes: 1 addition & 1 deletion src/atm_conv.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
You should have received a copy of the GNU General Public License
along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2013-2023 Forschungszentrum Juelich GmbH
Copyright (C) 2013-2025 Forschungszentrum Juelich GmbH
*/

/*!
Expand Down
75 changes: 41 additions & 34 deletions src/atm_dist.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
You should have received a copy of the GNU General Public License
along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2013-2024 Forschungszentrum Juelich GmbH
Copyright (C) 2013-2025 Forschungszentrum Juelich GmbH
*/

/*!
Expand All @@ -39,7 +39,7 @@ int main(
rqtdm[NQ], rvtdm, t0 =
0, x0[3], x1[3], x2[3], z1, *z1_old, z2, *z2_old, *work;

int f, init = 0, ip, iq, np;
int init = 0, np;

/* Allocate... */
ALLOC(atm1, atm_t, 1);
Expand Down Expand Up @@ -86,14 +86,21 @@ int main(

/* Read control parameters... */
read_ctl(argv[1], argc, argv, &ctl);
int ens = (int) scan_ctl(argv[1], argc, argv, "DIST_ENS", -1, "-999", NULL);
double p0 = P(scan_ctl(argv[1], argc, argv, "DIST_Z0", -1, "-1000", NULL));
double p1 = P(scan_ctl(argv[1], argc, argv, "DIST_Z1", -1, "1000", NULL));
double lat0 = scan_ctl(argv[1], argc, argv, "DIST_LAT0", -1, "-1000", NULL);
double lat1 = scan_ctl(argv[1], argc, argv, "DIST_LAT1", -1, "1000", NULL);
double lon0 = scan_ctl(argv[1], argc, argv, "DIST_LON0", -1, "-1000", NULL);
double lon1 = scan_ctl(argv[1], argc, argv, "DIST_LON1", -1, "1000", NULL);
double zscore =
const int ens =
(int) scan_ctl(argv[1], argc, argv, "DIST_ENS", -1, "-999", NULL);
const double p0 =
P(scan_ctl(argv[1], argc, argv, "DIST_Z0", -1, "-1000", NULL));
const double p1 =
P(scan_ctl(argv[1], argc, argv, "DIST_Z1", -1, "1000", NULL));
const double lat0 =
scan_ctl(argv[1], argc, argv, "DIST_LAT0", -1, "-1000", NULL);
const double lat1 =
scan_ctl(argv[1], argc, argv, "DIST_LAT1", -1, "1000", NULL);
const double lon0 =
scan_ctl(argv[1], argc, argv, "DIST_LON0", -1, "-1000", NULL);
const double lon1 =
scan_ctl(argv[1], argc, argv, "DIST_LON1", -1, "1000", NULL);
const double zscore =
scan_ctl(argv[1], argc, argv, "DIST_ZSCORE", -1, "-999", NULL);

/* Write info... */
Expand All @@ -112,7 +119,7 @@ int main(
"# $5 = absolute vertical distance (%s) [km]\n"
"# $6 = relative vertical distance (%s) [%%]\n",
argv[3], argv[3], argv[3], argv[3]);
for (iq = 0; iq < ctl.nq; iq++)
for (int iq = 0; iq < ctl.nq; iq++)
fprintf(out,
"# $%d = %s absolute difference (%s) [%s]\n"
"# $%d = %s relative difference (%s) [%%]\n",
Expand All @@ -121,7 +128,7 @@ int main(
fprintf(out, "# $%d = number of particles\n\n", 7 + 2 * ctl.nq);

/* Loop over file pairs... */
for (f = 4; f < argc; f += 2) {
for (int f = 4; f < argc; f += 2) {

/* Read atmopheric data... */
if (!read_atm(argv[f], &ctl, atm1) || !read_atm(argv[f + 1], &ctl, atm2))
Expand All @@ -132,7 +139,7 @@ int main(
ERRMSG("Different numbers of particles!");

/* Get time from filename... */
double t = time_from_filename(argv[f], ctl.atm_type < 2 ? 20 : 19);
const double t = time_from_filename(argv[f], ctl.atm_type < 2 ? 20 : 19);

/* Save initial time... */
if (!init) {
Expand All @@ -142,14 +149,14 @@ int main(

/* Init... */
np = 0;
for (ip = 0; ip < atm1->np; ip++) {
for (int ip = 0; ip < atm1->np; ip++) {
ahtd[ip] = avtd[ip] = rhtd[ip] = rvtd[ip] = 0;
for (iq = 0; iq < ctl.nq; iq++)
for (int iq = 0; iq < ctl.nq; iq++)
aqtd[iq * NP + ip] = rqtd[iq * NP + ip] = 0;
}

/* Loop over air parcels... */
for (ip = 0; ip < atm1->np; ip++) {
for (int ip = 0; ip < atm1->np; ip++) {

/* Check air parcel index... */
if (ctl.qnt_idx > 0
Expand Down Expand Up @@ -185,7 +192,7 @@ int main(
/* Calculate absolute transport deviations... */
ahtd[np] = DIST(x1, x2);
avtd[np] = z1 - z2;
for (iq = 0; iq < ctl.nq; iq++)
for (int iq = 0; iq < ctl.nq; iq++)
aqtd[iq * NP + np] = atm1->q[iq][ip] - atm2->q[iq][ip];

/* Calculate relative transport deviations... */
Expand All @@ -208,7 +215,7 @@ int main(
}

/* Get relative transport deviations... */
for (iq = 0; iq < ctl.nq; iq++)
for (int iq = 0; iq < ctl.nq; iq++)
rqtd[iq * NP + np] = 200. * (atm1->q[iq][ip] - atm2->q[iq][ip])
/ (fabs(atm1->q[iq][ip]) + fabs(atm2->q[iq][ip]));

Expand All @@ -229,11 +236,11 @@ int main(
if (zscore > 0 && np > 1) {

/* Get means and standard deviations of transport deviations... */
size_t n = (size_t) np;
double muh = gsl_stats_mean(ahtd, 1, n);
double muv = gsl_stats_mean(avtd, 1, n);
double sigh = gsl_stats_sd(ahtd, 1, n);
double sigv = gsl_stats_sd(avtd, 1, n);
const size_t n = (size_t) np;
const double muh = gsl_stats_mean(ahtd, 1, n);
const double muv = gsl_stats_mean(avtd, 1, n);
const double sigh = gsl_stats_sd(ahtd, 1, n);
const double sigv = gsl_stats_sd(avtd, 1, n);

/* Filter data... */
np = 0;
Expand All @@ -244,7 +251,7 @@ int main(
rhtd[np] = rhtd[i];
avtd[np] = avtd[i];
rvtd[np] = rvtd[i];
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtd[iq * NP + np] = aqtd[iq * NP + (int) i];
rqtd[iq * NP + np] = rqtd[iq * NP + (int) i];
}
Expand All @@ -258,7 +265,7 @@ int main(
rhtdm = gsl_stats_mean(rhtd, 1, (size_t) np);
avtdm = gsl_stats_mean(avtd, 1, (size_t) np);
rvtdm = gsl_stats_mean(rvtd, 1, (size_t) np);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_mean(&aqtd[iq * NP], 1, (size_t) np);
rqtdm[iq] = gsl_stats_mean(&rqtd[iq * NP], 1, (size_t) np);
}
Expand All @@ -267,7 +274,7 @@ int main(
rhtdm = gsl_stats_sd(rhtd, 1, (size_t) np);
avtdm = gsl_stats_sd(avtd, 1, (size_t) np);
rvtdm = gsl_stats_sd(rvtd, 1, (size_t) np);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_sd(&aqtd[iq * NP], 1, (size_t) np);
rqtdm[iq] = gsl_stats_sd(&rqtd[iq * NP], 1, (size_t) np);
}
Expand All @@ -276,7 +283,7 @@ int main(
rhtdm = gsl_stats_min(rhtd, 1, (size_t) np);
avtdm = gsl_stats_min(avtd, 1, (size_t) np);
rvtdm = gsl_stats_min(rvtd, 1, (size_t) np);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_min(&aqtd[iq * NP], 1, (size_t) np);
rqtdm[iq] = gsl_stats_min(&rqtd[iq * NP], 1, (size_t) np);
}
Expand All @@ -285,7 +292,7 @@ int main(
rhtdm = gsl_stats_max(rhtd, 1, (size_t) np);
avtdm = gsl_stats_max(avtd, 1, (size_t) np);
rvtdm = gsl_stats_max(rvtd, 1, (size_t) np);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_max(&aqtd[iq * NP], 1, (size_t) np);
rqtdm[iq] = gsl_stats_max(&rqtd[iq * NP], 1, (size_t) np);
}
Expand All @@ -294,7 +301,7 @@ int main(
rhtdm = gsl_stats_skew(rhtd, 1, (size_t) np);
avtdm = gsl_stats_skew(avtd, 1, (size_t) np);
rvtdm = gsl_stats_skew(rvtd, 1, (size_t) np);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_skew(&aqtd[iq * NP], 1, (size_t) np);
rqtdm[iq] = gsl_stats_skew(&rqtd[iq * NP], 1, (size_t) np);
}
Expand All @@ -303,7 +310,7 @@ int main(
rhtdm = gsl_stats_kurtosis(rhtd, 1, (size_t) np);
avtdm = gsl_stats_kurtosis(avtd, 1, (size_t) np);
rvtdm = gsl_stats_kurtosis(rvtd, 1, (size_t) np);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_kurtosis(&aqtd[iq * NP], 1, (size_t) np);
rqtdm[iq] = gsl_stats_kurtosis(&rqtd[iq * NP], 1, (size_t) np);
}
Expand All @@ -312,7 +319,7 @@ int main(
rhtdm = gsl_stats_absdev_m(rhtd, 1, (size_t) np, 0.0);
avtdm = gsl_stats_absdev_m(avtd, 1, (size_t) np, 0.0);
rvtdm = gsl_stats_absdev_m(rvtd, 1, (size_t) np, 0.0);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_absdev_m(&aqtd[iq * NP], 1, (size_t) np, 0.0);
rqtdm[iq] = gsl_stats_absdev_m(&rqtd[iq * NP], 1, (size_t) np, 0.0);
}
Expand All @@ -321,7 +328,7 @@ int main(
rhtdm = gsl_stats_median(rhtd, 1, (size_t) np);
avtdm = gsl_stats_median(avtd, 1, (size_t) np);
rvtdm = gsl_stats_median(rvtd, 1, (size_t) np);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_median(&aqtd[iq * NP], 1, (size_t) np);
rqtdm[iq] = gsl_stats_median(&rqtd[iq * NP], 1, (size_t) np);
}
Expand All @@ -330,7 +337,7 @@ int main(
rhtdm = gsl_stats_mad0(rhtd, 1, (size_t) np, work);
avtdm = gsl_stats_mad0(avtd, 1, (size_t) np, work);
rvtdm = gsl_stats_mad0(rvtd, 1, (size_t) np, work);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
aqtdm[iq] = gsl_stats_mad0(&aqtd[iq * NP], 1, (size_t) np, work);
rqtdm[iq] = gsl_stats_mad0(&rqtd[iq * NP], 1, (size_t) np, work);
}
Expand All @@ -340,7 +347,7 @@ int main(
/* Write output... */
fprintf(out, "%.2f %.2f %g %g %g %g", t, t - t0,
ahtdm, rhtdm, avtdm, rvtdm);
for (iq = 0; iq < ctl.nq; iq++) {
for (int iq = 0; iq < ctl.nq; iq++) {
fprintf(out, " ");
fprintf(out, ctl.qnt_format[iq], aqtdm[iq]);
fprintf(out, " ");
Expand Down
72 changes: 41 additions & 31 deletions src/atm_init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
You should have received a copy of the GNU General Public License
along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2013-2024 Forschungszentrum Juelich GmbH
Copyright (C) 2013-2025 Forschungszentrum Juelich GmbH
*/

/*!
Expand All @@ -32,8 +32,6 @@ int main(

ctl_t ctl;

gsl_rng *rng;

/* Allocate... */
ALLOC(atm, atm_t, 1);

Expand All @@ -43,38 +41,49 @@ int main(

/* Read control parameters... */
read_ctl(argv[1], argc, argv, &ctl);
double t0 = scan_ctl(argv[1], argc, argv, "INIT_T0", -1, "0", NULL);
double t1 = scan_ctl(argv[1], argc, argv, "INIT_T1", -1, "0", NULL);
double dt = scan_ctl(argv[1], argc, argv, "INIT_DT", -1, "1", NULL);
double z0 = scan_ctl(argv[1], argc, argv, "INIT_Z0", -1, "0", NULL);
double z1 = scan_ctl(argv[1], argc, argv, "INIT_Z1", -1, "0", NULL);
double dz = scan_ctl(argv[1], argc, argv, "INIT_DZ", -1, "1", NULL);
double lon0 = scan_ctl(argv[1], argc, argv, "INIT_LON0", -1, "0", NULL);
double lon1 = scan_ctl(argv[1], argc, argv, "INIT_LON1", -1, "0", NULL);
double dlon = scan_ctl(argv[1], argc, argv, "INIT_DLON", -1, "1", NULL);
double lat0 = scan_ctl(argv[1], argc, argv, "INIT_LAT0", -1, "0", NULL);
double lat1 = scan_ctl(argv[1], argc, argv, "INIT_LAT1", -1, "0", NULL);
double dlat = scan_ctl(argv[1], argc, argv, "INIT_DLAT", -1, "1", NULL);
double st = scan_ctl(argv[1], argc, argv, "INIT_ST", -1, "0", NULL);
double sz = scan_ctl(argv[1], argc, argv, "INIT_SZ", -1, "0", NULL);
double slon = scan_ctl(argv[1], argc, argv, "INIT_SLON", -1, "0", NULL);
double slat = scan_ctl(argv[1], argc, argv, "INIT_SLAT", -1, "0", NULL);
double sx = scan_ctl(argv[1], argc, argv, "INIT_SX", -1, "0", NULL);
double ut = scan_ctl(argv[1], argc, argv, "INIT_UT", -1, "0", NULL);
double uz = scan_ctl(argv[1], argc, argv, "INIT_UZ", -1, "0", NULL);
double ulon = scan_ctl(argv[1], argc, argv, "INIT_ULON", -1, "0", NULL);
double ulat = scan_ctl(argv[1], argc, argv, "INIT_ULAT", -1, "0", NULL);
int even =
const double t0 = scan_ctl(argv[1], argc, argv, "INIT_T0", -1, "0", NULL);
const double t1 = scan_ctl(argv[1], argc, argv, "INIT_T1", -1, "0", NULL);
const double dt = scan_ctl(argv[1], argc, argv, "INIT_DT", -1, "1", NULL);
const double z0 = scan_ctl(argv[1], argc, argv, "INIT_Z0", -1, "0", NULL);
const double z1 = scan_ctl(argv[1], argc, argv, "INIT_Z1", -1, "0", NULL);
const double dz = scan_ctl(argv[1], argc, argv, "INIT_DZ", -1, "1", NULL);
const double lon0 =
scan_ctl(argv[1], argc, argv, "INIT_LON0", -1, "0", NULL);
const double lon1 =
scan_ctl(argv[1], argc, argv, "INIT_LON1", -1, "0", NULL);
const double dlon =
scan_ctl(argv[1], argc, argv, "INIT_DLON", -1, "1", NULL);
const double lat0 =
scan_ctl(argv[1], argc, argv, "INIT_LAT0", -1, "0", NULL);
const double lat1 =
scan_ctl(argv[1], argc, argv, "INIT_LAT1", -1, "0", NULL);
const double dlat =
scan_ctl(argv[1], argc, argv, "INIT_DLAT", -1, "1", NULL);
const double st = scan_ctl(argv[1], argc, argv, "INIT_ST", -1, "0", NULL);
const double sz = scan_ctl(argv[1], argc, argv, "INIT_SZ", -1, "0", NULL);
const double slon =
scan_ctl(argv[1], argc, argv, "INIT_SLON", -1, "0", NULL);
const double slat =
scan_ctl(argv[1], argc, argv, "INIT_SLAT", -1, "0", NULL);
const double sx = scan_ctl(argv[1], argc, argv, "INIT_SX", -1, "0", NULL);
const double ut = scan_ctl(argv[1], argc, argv, "INIT_UT", -1, "0", NULL);
const double uz = scan_ctl(argv[1], argc, argv, "INIT_UZ", -1, "0", NULL);
const double ulon =
scan_ctl(argv[1], argc, argv, "INIT_ULON", -1, "0", NULL);
const double ulat =
scan_ctl(argv[1], argc, argv, "INIT_ULAT", -1, "0", NULL);
const int even =
(int) scan_ctl(argv[1], argc, argv, "INIT_EVENLY", -1, "0", NULL);
int rep = (int) scan_ctl(argv[1], argc, argv, "INIT_REP", -1, "1", NULL);
double m = scan_ctl(argv[1], argc, argv, "INIT_MASS", -1, "0", NULL);
double vmr = scan_ctl(argv[1], argc, argv, "INIT_VMR", -1, "0", NULL);
double bellrad =
const int rep =
(int) scan_ctl(argv[1], argc, argv, "INIT_REP", -1, "1", NULL);
const double m = scan_ctl(argv[1], argc, argv, "INIT_MASS", -1, "0", NULL);
const double vmr = scan_ctl(argv[1], argc, argv, "INIT_VMR", -1, "0", NULL);
const double bellrad =
scan_ctl(argv[1], argc, argv, "INIT_BELLRAD", -1, "0", NULL);

/* Initialize random number generator... */
gsl_rng_env_setup();
rng = gsl_rng_alloc(gsl_rng_default);
gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

/* Create grid... */
for (double t = t0; t <= t1; t += dt)
Expand Down Expand Up @@ -111,7 +120,8 @@ int main(
double x0[3], x1[3];
geo2cart(0.0, 0.5 * (lon0 + lon1), 0.5 * (lat0 + lat1), x0);
geo2cart(0.0, atm->lon[atm->np], atm->lat[atm->np], x1);
double rad = RE * acos(DOTP(x0, x1) / NORM(x0) / NORM(x1));
const double rad =
RE * acos(DOTP(x0, x1) / NORM(x0) / NORM(x1));
if (rad > bellrad)
continue;
if (ctl.qnt_m >= 0)
Expand Down
Loading

0 comments on commit a25a378

Please sign in to comment.