AIRS Code Collection
Macros | Functions
rayt.c File Reference

2-D gravity wave ray-tracing code. More...

#include "libairs.h"

Go to the source code of this file.

Macros

#define NZ   1000
 

Functions

double buoyancy (double z0, double p0, double t0, double z1, double p1, double t1)
 
double scale_height (double t)
 
double temp2theta (double p, double t)
 
int main (int argc, char *argv[])
 

Detailed Description

2-D gravity wave ray-tracing code.

Definition in file rayt.c.

Macro Definition Documentation

◆ NZ

#define NZ   1000

Definition at line 33 of file rayt.c.

Function Documentation

◆ buoyancy()

double buoyancy ( double  z0,
double  p0,
double  t0,
double  z1,
double  p1,
double  t1 
)

Definition at line 231 of file rayt.c.

237 {
238
239 /* Get potential temperature... */
240 const double theta0 = temp2theta(p0, t0);
241 const double theta1 = temp2theta(p1, t1);
242
243 /* Get buoyancy frequency... */
244 return sqrt(G0 / (0.5 * (theta0 + theta1)) * (theta1 - theta0) /
245 ((z1 - z0) * 1e3));
246}
G0
#define G0
Standard gravity [m/s^2].
Definition: jurassic.h:274
temp2theta
double temp2theta(double p, double t)
Definition: rayt.c:258
Here is the call graph for this function:

◆ scale_height()

double scale_height ( double  t)

Definition at line 250 of file rayt.c.

251 {
252
253 return 29.26 * t / 1e3;
254}

◆ temp2theta()

double temp2theta ( double  p,
double  t 
)

Definition at line 258 of file rayt.c.

260 {
261
262 return t * pow(P0 / p, 0.286);
263}
P0
#define P0
Standard pressure [hPa].
Definition: jurassic.h:309

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 61 of file rayt.c.

63 {
64
65 FILE *in;
66
67 const double dzw = 5 * 1e3;
68
69 static double f0, z[NZ], u[NZ], urel[NZ], v[NZ], bf[NZ], bf2[NZ], H[NZ],
70 frel[NZ], osign[NZ], f1[NZ], f2[NZ], delta[NZ], a2[NZ], m[NZ], dxdz[NZ],
71 cgz[NZ], dz, path[NZ], tim[NZ], p[NZ], t[NZ], fgb;
72
73 static int izcrit, izrefl, nz;
74
75 /* Check arguments... */
76 if (argc != 8)
77 ERRMSG("Give parameters: <atm.tab> <z_launch> <mode> "
78 "<t_gb | lz_launch> <lx> <lat> <direct>");
79
80 /* Get launch level... */
81 const double z0 = atof(argv[2]);
82 const double lat = atof(argv[6]);
83 const double alpha = atof(argv[7]);
84
85 /* Read atmosphere above launch level... */
86 if (!(in = fopen(argv[1], "r")))
87 ERRMSG("Cannot open atmospheric data file!");
88 while (fscanf
89 (in, "%lg %lg %lg %lg %lg", &z[nz], &p[nz], &t[nz], &u[nz], &v[nz])
90 == 5)
91 if (z[nz] >= z0) {
92 u[nz] =
93 cos(alpha * M_PI / 180.) * u[nz] + sin(alpha * M_PI / 180.) * v[nz];
94 if ((++nz) > NZ)
95 ERRMSG("Too many altitude levels!");
96 }
97 fclose(in);
98
99 /* Compute scale height and buoyancy frequency... */
100 for (int iz = 0; iz < nz; iz++) {
101 if (iz < nz - 1)
102 bf[iz] = buoyancy(z[iz], p[iz], t[iz], z[iz + 1], p[iz + 1], t[iz + 1]);
103 else
104 bf[iz] = bf[iz - 1];
105 H[iz] = scale_height(t[iz]) * 1e3;
106 z[iz] *= 1e3;
107 }
108
109 /* Smooth N profile... */
110 for (int iz = 0; iz < nz; iz++) {
111 bf2[iz] = 0;
112 double wsum = 0;
113 for (int iz2 = 0; iz2 < nz; iz2++) {
114 if (!gsl_finite(bf[iz2]) ||
115 !gsl_finite(bf[GSL_MAX(iz2 - 1, 0)]) ||
116 !gsl_finite(bf[GSL_MIN(iz2 + 1, nz - 1)]))
117 continue;
118 const double w =
119 (fabs(z[iz] - z[iz2]) < dzw) ? 1.0 - fabs(z[iz] - z[iz2]) / dzw : 0.0;
120 bf2[iz] += w * bf[iz2];
121 wsum += w;
122 }
123 bf2[iz] /= wsum;
124 }
125 for (int iz = 0; iz < nz; iz++)
126 bf[iz] = bf2[iz];
127
128 /* Get horizontal wavenumber... */
129 const double k = 2 * M_PI / (atof(argv[5]) * 1e3);
130
131 /* Get minimum gravity wave frequency (Coriolis parameter)... */
132 const double omin = 2 * 2 * M_PI / 86400. * sin(lat / 180. * M_PI);
133
134 /* Get initial frequencies... */
135 if (argv[3][0] == 't') {
136
137 /* Get ground-based frequency... */
138 fgb = 2 * M_PI / (atof(argv[4]) * 60.);
139
140 /* Get intrinsic frequency at launch level... */
141 f0 = fgb - k * u[0];
142
143 } else if (argv[3][0] == 'l') {
144
145 /* Get vertical wavenumber... */
146 const double m0 = 2 * M_PI / (atof(argv[4]) * 1e3);
147
148 /* Get intrinsic frequency at launch level... */
149 f0 =
150 sqrt((bf[0] * bf[0] * k * k +
151 omin * omin * (m0 * m0 + 0.25 / (H[0] * H[0])))
152 / (m0 * m0 + k * k + 0.25 / (H[0] * H[0])));
153
154 /* Get ground-based frequency... */
155 fgb = f0 + k * u[0];
156
157 } else
158 ERRMSG("Set <mode> to 't_gb' or 'lz_launch'!");
159
160 /* Loop over layers... */
161 for (int iz = 0; iz < nz; iz++) {
162 urel[iz] = u[iz] - u[0];
163 frel[iz] = f0 - k * urel[iz];
164 osign[iz] = frel[iz] / fabs(frel[iz]);
165 f1[iz] = (bf[iz] * bf[iz] - frel[iz] * frel[iz]) / frel[iz];
166 f2[iz] = (frel[iz] * frel[iz] - omin * omin) / frel[iz];
167 delta[iz] = k * k * (1 + f1[iz] / f2[iz]);
168 a2[iz] = 1. / 4. / (H[iz] * H[iz]);
169 m[iz] = (-osign[iz]) * k * sqrt((f1[iz] / f2[iz]) - (a2[iz] / (k * k)));
170 dxdz[iz] = (u[iz] * delta[iz] + k * f1[iz]) / (-1 * m[iz] * f2[iz]);
171 dz = z[1] - z[0];
172 cgz[iz] = f2[iz] * (-1. * m[iz]) / (k * k + m[iz] * m[iz] + a2[iz]);
173 }
174
175 /* Integrate via trapezoidal rule... */
176 for (int iz = 1; iz < nz; iz++) {
177 path[iz] = path[iz - 1] + dz * .5 * (dxdz[iz - 1] + dxdz[iz]);
178 tim[iz] = tim[iz - 1] + dz * 2. / (cgz[iz - 1] + cgz[iz]);
179 }
180
181 /* Find critical level... */
182 for (izcrit = 0; izcrit < nz; izcrit++)
183 if (f0 / fabs(f0) * frel[izcrit] / fabs(omin) <= 1)
184 break;
185
186 /* Find trapping/reflection level... */
187 for (izrefl = 0; izrefl < nz; izrefl++) {
188 const double costh = fabs(f0 - k * urel[izrefl])
189 / sqrt(bf[izrefl] * bf[izrefl]
190 * (1 -
191 (1 -
192 (omin / bf[izrefl]) * (omin / bf[izrefl])) / (k * k /
193 a2[izrefl] +
194 1)));
195 if (costh >= 1.0)
196 break;
197 }
198
199 /* Filter data... */
200 for (int iz = 0; iz < nz; iz++)
201 if (iz >= izcrit || iz >= izrefl)
202 path[iz] = tim[iz] = m[iz] = frel[iz] = cgz[iz] = sqrt(-1.0);
203
204 /* Write output... */
205 printf("# $1 = latitude [deg]\n"
206 "# $2 = altitude [km]\n"
207 "# $3 = pressure [hPa]\n"
208 "# $4 = temperature [K]\n"
209 "# $5 = potential temperature [K]\n"
210 "# $6 = wind speed [m/s]\n"
211 "# $7 = buoyancy frequency [1/s]\n"
212 "# $8 = scale height [km]\n"
213 "# $9 = horizontal distance [km]\n"
214 "# $10 = propagation time [min]\n"
215 "# $11 = vertical wavelength [km]\n"
216 "# $12 = wave period [min]\n"
217 "# $13 = vertical group velocity [m/s]\n\n");
218 for (int iz = 0; iz < nz; iz++)
219 printf("%g %g %g %g %g %g %g %g %g %g %g %g %g\n",
220 lat, z[iz] / 1e3, p[iz], t[iz], temp2theta(p[iz], t[iz]), u[iz],
221 bf[iz], H[iz] / 1e3, path[iz] / 1e3, tim[iz] / 60,
222 fabs(2 * M_PI / m[iz] / 1e3), 2. * M_PI / frel[iz] / 60., cgz[iz]);
223 printf("\n# z_crit= %g km\n# z_refl= %g km\n",
224 z[izcrit - 1] / 1e3, z[izrefl - 1] / 1e3);
225
226 return EXIT_SUCCESS;
227}
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
scale_height
double scale_height(double t)
Definition: rayt.c:250
NZ
#define NZ
Definition: rayt.c:33
buoyancy
double buoyancy(double z0, double p0, double t0, double z1, double p1, double t1)
Definition: rayt.c:231
Here is the call graph for this function:
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4