41 {
42
44
46
48
50
51
53
54
55 if (argc < 4)
56 ERRMSG(
"Missing or invalid command-line arguments.\n\n"
57 "Usage: met_check_dt <ctl> <dt_file> <met> [KEY VALUE ...]\n\n"
58 "Use -h for full help.");
59
60
64
65
67 const double kx =
scan_ctl(argv[1], argc, argv,
"KX", -1,
"50.0", NULL);
68 const double kz =
scan_ctl(argv[1], argc, argv,
"KZ", -1,
"0.1", NULL);
69 const double dx = 1e3 *
scan_ctl(argv[1], argc, argv,
"DX", -1,
"", NULL);
70 const double c_max =
scan_ctl(argv[1], argc, argv,
"CMAX", -1,
"0.5", NULL);
71 const double n_max =
scan_ctl(argv[1], argc, argv,
"NMAX", -1,
"0.3", NULL);
72
73
75
76
78 ERRMSG(
"Cannot open file!");
79
80
81 FILE *out;
82 LOG(1,
"Write time step data file: %s", argv[2]);
83 if (!(out = fopen(argv[2], "w")))
84 ERRMSG(
"Cannot create file!");
85
86
87 fprintf(out,
88 "# $1 = height [km]\n"
89 "# $2 = time step for horizontal advection [s]\n"
90 "# $3 = time step for vertical advection [s]\n"
91 "# $4 = time step for horizontal diffusion [s]\n"
92 "# $5 = time step for vertical diffusion [s]\n\n");
93
94
95 for (
int ip = 1; ip < met->
np - 1; ip++) {
96
97
98 double dt_x_min = 1e100;
99 double dt_p_min = 1e100;
100 double dt_dx_min = 1e100;
101 double dt_dp_min = 1e100;
102
103
104#pragma omp parallel for default(shared) collapse(2) reduction(min:dt_x_min,dt_p_min,dt_dx_min,dt_dp_min)
105 for (
int ix = 0; ix < met->
nx; ix++)
106 for (
int iy = 1; iy < met->
ny - 1; iy++) {
107
108
109 const double vh =
110 sqrt(
SQR(met->
u[ix][iy][ip]) +
SQR(met->
v[ix][iy][ip]));
111 const double dt_x = fabs(c_max * dx / vh);
112 if (vh != 0)
113 dt_x_min =
MIN(dt_x, dt_x_min);
114
115 const double dp = 0.5 * fabs(met->
p[ip + 1] - met->
p[ip - 1]);
116 const double dt_p = fabs(c_max * dp / met->
w[ix][iy][ip]);
117 if (met->
w[ix][iy][ip] != 0)
118 dt_p_min =
MIN(dt_p, dt_p_min);
119
120
121 const double dt_dx = 0.5 *
SQR(n_max * dx) / kx;
122 dt_dx_min =
MIN(dt_dx, dt_dx_min);
123
124 const double dt_dp =
125 0.5 *
SQR(n_max * dp) / (
SQR(met->
p[ip] / (100. *
H0)) * kz);
126 dt_dp_min =
MIN(dt_dp, dt_dp_min);
127 }
128
129
130 fprintf(out,
"%g %g %g %g %g\n",
Z(met->
p[ip]), dt_x_min, dt_p_min,
131 dt_dx_min, dt_dp_min);
132 }
133
134
135 fclose(out);
136
137
138 free(clim);
139 free(met);
140 free(dd);
141
142 return EXIT_SUCCESS;
143}
double scan_ctl(const char *filename, int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Scans a control file or command-line arguments for a specified variable.
void mptrac_read_clim(const ctl_t *ctl, clim_t *clim)
Reads various climatological data and populates the given climatology structure.
int mptrac_read_met(const char *filename, const ctl_t *ctl, const clim_t *clim, met_t *met, dd_t *dd)
Reads meteorological data from a file, supporting multiple formats and MPI broadcasting.
void mptrac_read_ctl(const char *filename, int argc, char *argv[], ctl_t *ctl)
Reads control parameters from a configuration file and populates the given structure.
#define H0
Scale height [km].
#define MIN(a, b)
Macro to determine the minimum of two values.
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
#define USAGE
Print usage information on -h or --help.
#define Z(p)
Convert pressure to altitude.
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
#define SQR(x)
Compute the square of a value.
#define LOG(level,...)
Print a log message with a specified logging level.
Domain decomposition data structure.
float w[EX][EY][EP]
Vertical velocity [hPa/s].
int nx
Number of longitudes.
int ny
Number of latitudes.
int np
Number of pressure levels.
float u[EX][EY][EP]
Zonal wind [m/s].
float v[EX][EY][EP]
Meridional wind [m/s].
double p[EP]
Pressure levels [hPa].