1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7 #include <cmath>
8
9 #ifdef HAVE_CONFIG_H
10 #include "config.h"
11 #endif
12
13 #include "constants.h"
14
15 #define SQUARE_RADIUS (-PlanetRadius * PlanetRadius)
16
17 void
dv2ps(const double * restrict div,double * restrict pot,long nlev,long ntr)18 dv2ps(const double *restrict div, double *restrict pot, long nlev, long ntr)
19 {
20 long l, m, n;
21 double fact;
22
23 for (l = 0; l < nlev; l++)
24 {
25 /* m == 0 */
26 *pot++ = 0.0;
27 *pot++ = 0.0;
28 div += 2;
29
30 for (n = 1; n <= ntr; n++)
31 {
32 fact = SQUARE_RADIUS / (n * n + n);
33 *pot++ = *div++ * fact;
34 *pot++ = *div++ * fact;
35 }
36
37 /* m >= 0 */
38 for (m = 1; m <= ntr; m++)
39 for (n = m; n <= ntr; n++)
40 {
41 fact = SQUARE_RADIUS / (n * n + n);
42 *pot++ = *div++ * fact;
43 *pot++ = *div++ * fact;
44 }
45 }
46 }
47
48 static void
dv2uv_kernel(const double * d,const double * o,double * restrict u,double * restrict v,const double * f,const double * g,long nt)49 dv2uv_kernel(const double *d, const double *o, double *restrict u, double *restrict v, const double *f, const double *g, long nt)
50 {
51 // d(nsp,nlev), o(nsp,nlev) ! divergence, vorticity
52 // u(nsp,nlev), v(nsp,nlev) ! zonal wind, meridional wind
53 // f(nsp/2) , g(nsp/2) ! factor tables
54
55 long i = 0;
56 for (long m = 0; m < nt - 1; m++)
57 {
58 // n = m
59 if (m == 0)
60 {
61 *u++ = -g[i + 1] * o[2 * (i + 1)];
62 *u++ = -g[i + 1] * o[2 * (i + 1) + 1];
63 *v++ = g[i + 1] * d[2 * (i + 1)];
64 *v++ = g[i + 1] * d[2 * (i + 1) + 1];
65 }
66 else
67 {
68 *u++ = -f[i] * d[2 * i + 1] - g[i + 1] * o[2 * (i + 1)];
69 *u++ = f[i] * d[2 * i] - g[i + 1] * o[2 * (i + 1) + 1];
70 *v++ = -f[i] * o[2 * i + 1] + g[i + 1] * d[2 * (i + 1)];
71 *v++ = f[i] * o[2 * i] + g[i + 1] * d[2 * (i + 1) + 1];
72 }
73 ++i;
74
75 // m < n < nt-1
76 for (long n = m + 1; n < nt - 1; n++)
77 {
78 *u++ = g[i] * o[2 * (i - 1)] - f[i] * d[2 * i + 1] - g[i + 1] * o[2 * (i + 1)];
79 *u++ = g[i] * o[2 * (i - 1) + 1] + f[i] * d[2 * i] - g[i + 1] * o[2 * (i + 1) + 1];
80 *v++ = -g[i] * d[2 * (i - 1)] - f[i] * o[2 * i + 1] + g[i + 1] * d[2 * (i + 1)];
81 *v++ = -g[i] * d[2 * (i - 1) + 1] + f[i] * o[2 * i] + g[i + 1] * d[2 * (i + 1) + 1];
82 ++i;
83 }
84
85 // n = nt-1
86 *u++ = g[i] * o[2 * (i - 1)] - f[i] * d[2 * i + 1];
87 *u++ = g[i] * o[2 * (i - 1) + 1] + f[i] * d[2 * i];
88 *v++ = -g[i] * d[2 * (i - 1)] - f[i] * o[2 * i + 1];
89 *v++ = -g[i] * d[2 * (i - 1) + 1] + f[i] * o[2 * i];
90 ++i;
91
92 // n = nt
93 *u++ = g[i] * o[2 * (i - 1)];
94 *u++ = g[i] * o[2 * (i - 1) + 1];
95 *v++ = -g[i] * d[2 * (i - 1)];
96 *v++ = -g[i] * d[2 * (i - 1) + 1];
97 ++i;
98 }
99
100 // m = nt-1 and n = nt-1
101 *u++ = -f[i] * d[2 * i + 1];
102 *u++ = f[i] * d[2 * i];
103 *v++ = -f[i] * o[2 * i + 1];
104 *v++ = f[i] * o[2 * i];
105 ++i;
106
107 // m = nt-1 and n = nt
108 *u++ = g[i] * o[2 * (i - 1)];
109 *u++ = g[i] * o[2 * (i - 1) + 1];
110 *v++ = -g[i] * d[2 * (i - 1)];
111 *v++ = -g[i] * d[2 * (i - 1) + 1];
112 ++i;
113
114 // m = nt and n = nt
115 *u++ = 0.0;
116 *u++ = 0.0;
117 *v++ = 0.0;
118 *v++ = 0.0;
119 }
120
121 void
dv2uv(const double * d,const double * o,double * u,double * v,const double * f,const double * g,long nt,long nsp,long nlev)122 dv2uv(const double *d, const double *o, double *u, double *v, const double *f, const double *g, long nt, long nsp, long nlev)
123 {
124 // d(nsp,nlev), o(nsp,nlev) ! divergence, vorticity
125 // u(nsp,nlev), v(nsp,nlev) ! zonal wind, meridional wind
126 // f(nsp/2) , g(nsp/2) ! factor tables
127
128 #if defined(_OPENMP)
129 #pragma omp parallel for default(none) shared(d, o, u, v, f, g, nt, nsp, nlev)
130 #endif
131 for (long l = 0; l < nlev; l++)
132 {
133 dv2uv_kernel(d + l * nsp, o + l * nsp, u + l * nsp, v + l * nsp, f, g, nt);
134 }
135 }
136
137 /*
138 void scaluv(double *fu, const double *rclat, long nlat, long lot)
139 {
140 long l, lat;
141 double *ful;
142
143 #if defined (_OPENMP)
144 #pragma omp parallel for default(shared) private(lat, ful)
145 #endif
146 for ( l = 0; l < lot; l++ )
147 {
148 ful = fu + l*nlat;
149 for ( lat = 0; lat < nlat; lat++ )
150 {
151 ful[lat] = rclat[lat];
152 }
153 }
154 }
155 */
156
157 void
scaluv(double * fu,const double * rclat,long nlat,long lot)158 scaluv(double *fu, const double *rclat, long nlat, long lot)
159 {
160 for (long l = 0; l < lot; l++)
161 for (long lat = 0; lat < nlat; lat++)
162 {
163 *fu *= rclat[lat];
164 fu++;
165 }
166 }
167
168 void
uv2dv(const double * fu,const double * fv,double * sd,double * sv,const double * pol2,const double * pol3,long klev,long nlat,long nt)169 uv2dv(const double *fu, const double *fv, double *sd, double *sv, const double *pol2, const double *pol3, long klev, long nlat,
170 long nt)
171 {
172 long nsp2 = (nt + 1) * (nt + 2);
173 long nfc = (nt + 1) * 2;
174
175 #if defined(_OPENMP)
176 #pragma omp parallel for default(shared)
177 #endif
178 for (long lev = 0; lev < klev; lev++)
179 {
180 auto po2 = pol2;
181 auto po3 = pol3;
182 auto ful = fu + lev * nfc * nlat;
183 auto fvl = fv + lev * nfc * nlat;
184 auto sdl = sd + lev * nsp2;
185 auto svl = sv + lev * nsp2;
186 for (long jmm = 0; jmm <= nt; jmm++)
187 {
188 for (long jfc = jmm; jfc <= nt; jfc++)
189 {
190 auto ufr = ful;
191 auto ufi = ful + nlat;
192 auto vfr = fvl;
193 auto vfi = fvl + nlat;
194 double dir = 0.0;
195 double dii = 0.0;
196 double vor = 0.0;
197 double voi = 0.0;
198 for (long lat = 0; lat < nlat; lat++)
199 {
200 dir += vfr[lat] * po2[lat] - ufi[lat] * po3[lat];
201 dii += vfi[lat] * po2[lat] + ufr[lat] * po3[lat];
202 vor -= ufr[lat] * po2[lat] + vfi[lat] * po3[lat];
203 voi -= ufi[lat] * po2[lat] - vfr[lat] * po3[lat];
204 }
205 *sdl++ = dir;
206 *sdl++ = dii;
207 *svl++ = vor;
208 *svl++ = voi;
209 po2 += nlat;
210 po3 += nlat;
211 }
212 ful += 2 * nlat;
213 fvl += 2 * nlat;
214 }
215 }
216 }
217
218 void
geninx(long ntr,double * f,double * g)219 geninx(long ntr, double *f, double *g)
220 {
221 for (long m = 0; m <= ntr; m++)
222 {
223 long m2 = m * m;
224 for (long n = m; n <= ntr; n++)
225 {
226 long n2 = n * n;
227 if (n)
228 {
229 *g++ = -PlanetRadius / n * std::sqrt((double) (n2 - m2) / (double) (4 * n2 - 1));
230 *f++ = -PlanetRadius * m / (double) (n2 + n);
231 }
232 else
233 {
234 *g++ = 0.0;
235 *f++ = 0.0;
236 }
237 }
238 }
239 }
240