1 | ! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt |
---|
2 | ! SPDX-License-Identifier: GPL-3.0-or-later |
---|
3 | |
---|
4 | ! Changes to the routines by A. Stohl |
---|
5 | ! xi,xi0,eta,eta0 are double precision variables to avoid problems |
---|
6 | ! at poles |
---|
7 | |
---|
8 | module cmapf_mod |
---|
9 | |
---|
10 | use par_mod, only: dp |
---|
11 | |
---|
12 | implicit none |
---|
13 | private |
---|
14 | |
---|
15 | public :: cc2gll, cll2xy, cgszll, cxy2ll, stlmbr, stcm2p |
---|
16 | |
---|
17 | real,parameter :: rearth=6371.2, almst1=.9999999 |
---|
18 | |
---|
19 | real,parameter :: pi=3.14159265358979 |
---|
20 | real,parameter :: radpdg=pi/180., dgprad=180./pi |
---|
21 | |
---|
22 | contains |
---|
23 | |
---|
24 | subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg) |
---|
25 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
26 | |
---|
27 | use par_mod, only: dp |
---|
28 | |
---|
29 | implicit none |
---|
30 | |
---|
31 | real :: strcmp(9), xlat, xlong, ue, vn, ug, vg |
---|
32 | |
---|
33 | real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot |
---|
34 | |
---|
35 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
36 | if (xlat.gt.89.985) then |
---|
37 | !* North polar meteorological orientation: "north" along prime meridian |
---|
38 | rot = - strcmp(1) * along + xlong - 180. |
---|
39 | elseif (xlat.lt.-89.985) then |
---|
40 | !* South polar meteorological orientation: "north" along prime meridian |
---|
41 | rot = - strcmp(1) * along - xlong |
---|
42 | else |
---|
43 | rot = - strcmp(1) * along |
---|
44 | endif |
---|
45 | slong = sin( radpdg * rot ) |
---|
46 | clong = cos( radpdg * rot ) |
---|
47 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
48 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
49 | ug = ypolg * ue + xpolg * vn |
---|
50 | vg = ypolg * vn - xpolg * ue |
---|
51 | return |
---|
52 | end subroutine cc2gll |
---|
53 | |
---|
54 | subroutine ccrvll (strcmp, xlat,xlong, gx,gy) |
---|
55 | !* Written on 9/20/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
56 | |
---|
57 | use par_mod, only: dp |
---|
58 | |
---|
59 | implicit none |
---|
60 | |
---|
61 | real(kind=dp) :: xpolg,ypolg,temp,along,slong,clong,ctemp, curv |
---|
62 | real :: strcmp(9), xlat, xlong, gx, gy |
---|
63 | |
---|
64 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
65 | slong = sin( radpdg * strcmp(1) * along) |
---|
66 | clong = cos( radpdg * strcmp(1) * along) |
---|
67 | xpolg = - slong * strcmp(5) + clong * strcmp(6) |
---|
68 | ypolg = clong * strcmp(5) + slong * strcmp(6) |
---|
69 | temp = sin(radpdg * xlat) |
---|
70 | ctemp = cos(radpdg * xlat) |
---|
71 | curv = (strcmp(1) - temp) / ctemp / rearth |
---|
72 | gx = curv * xpolg |
---|
73 | gy = curv * ypolg |
---|
74 | return |
---|
75 | end subroutine ccrvll |
---|
76 | |
---|
77 | subroutine ccrvxy (strcmp, x,y, gx,gy) |
---|
78 | !* Written on 9/20/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
79 | |
---|
80 | use par_mod, only: dp |
---|
81 | |
---|
82 | implicit none |
---|
83 | |
---|
84 | real :: strcmp(9), x, y, gx, gy |
---|
85 | real(kind=dp) :: xpolg,ypolg,temp,ymerc,efact,curv |
---|
86 | |
---|
87 | temp = strcmp(1) * strcmp(7) /rearth |
---|
88 | xpolg = strcmp(6) + temp * (strcmp(3) - x) |
---|
89 | ypolg = strcmp(5) + temp * (strcmp(4) - y) |
---|
90 | temp = sqrt ( xpolg ** 2 + ypolg ** 2 ) |
---|
91 | if (temp.gt.0.) then |
---|
92 | ymerc = - log( temp) /strcmp(1) |
---|
93 | efact = exp(ymerc) |
---|
94 | curv = ( (strcmp(1) - 1.d0) * efact + & |
---|
95 | (strcmp(1) + 1.d0) / efact ) & |
---|
96 | * .5d0 / rearth |
---|
97 | gx = xpolg * curv / temp |
---|
98 | gy = ypolg * curv / temp |
---|
99 | else |
---|
100 | if (abs(strcmp(1)) .eq. 1.) then |
---|
101 | gx = 0. |
---|
102 | gy = 0. |
---|
103 | else |
---|
104 | gx = 1./rearth |
---|
105 | gy = 1./rearth |
---|
106 | endif |
---|
107 | endif |
---|
108 | return |
---|
109 | end subroutine ccrvxy |
---|
110 | |
---|
111 | subroutine cg2cll (strcmp, xlat,xlong, ug,vg, ue,vn) |
---|
112 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
113 | |
---|
114 | use par_mod, only: dp |
---|
115 | |
---|
116 | implicit none |
---|
117 | |
---|
118 | real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot |
---|
119 | real :: strcmp(9), xlat, xlong, ug, vg, ue, vn |
---|
120 | |
---|
121 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
122 | if (xlat.gt.89.985) then |
---|
123 | !* North polar meteorological orientation: "north" along prime meridian |
---|
124 | rot = - strcmp(1) * along + xlong - 180. |
---|
125 | elseif (xlat.lt.-89.985) then |
---|
126 | !* South polar meteorological orientation: "north" along prime meridian |
---|
127 | rot = - strcmp(1) * along - xlong |
---|
128 | else |
---|
129 | rot = - strcmp(1) * along |
---|
130 | endif |
---|
131 | slong = sin( radpdg * rot ) |
---|
132 | clong = cos( radpdg * rot ) |
---|
133 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
134 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
135 | ue = ypolg * ug - xpolg * vg |
---|
136 | vn = ypolg * vg + xpolg * ug |
---|
137 | return |
---|
138 | end subroutine cg2cll |
---|
139 | |
---|
140 | subroutine cg2cxy (strcmp, x,y, ug,vg, ue,vn) |
---|
141 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
142 | |
---|
143 | use par_mod, only: dp |
---|
144 | |
---|
145 | implicit none |
---|
146 | |
---|
147 | real :: strcmp(9) , x, y, ug, vg, ue, vn |
---|
148 | |
---|
149 | real :: clong, radial, rot, slong, xlat, xlong |
---|
150 | real(kind=dp) :: xpolg,ypolg,temp,xi0,eta0,xi,eta |
---|
151 | |
---|
152 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
153 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
154 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
155 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
156 | radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta) |
---|
157 | if (radial.gt.strcmp(8)) then |
---|
158 | !* Case north of 89 degrees. Meteorological wind direction definition |
---|
159 | !* changes. |
---|
160 | call cnxyll(strcmp, xi,eta, xlat,xlong) |
---|
161 | !* North polar meteorological orientation: "north" along prime meridian |
---|
162 | rot = strcmp(1) * (xlong - strcmp(2)) - xlong - 180. |
---|
163 | slong = - sin( radpdg * rot ) |
---|
164 | clong = cos( radpdg * rot ) |
---|
165 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
166 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
167 | else if (radial.lt.strcmp(9)) then |
---|
168 | !* Case south of -89 degrees. Meteorological wind direction definition |
---|
169 | !* changes. |
---|
170 | call cnxyll(strcmp, xi,eta, xlat,xlong) |
---|
171 | !* South polar meteorological orientation: "north" along prime meridian |
---|
172 | rot = strcmp(1) * (xlong - strcmp(2)) + xlong |
---|
173 | slong = - sin( radpdg * rot ) |
---|
174 | clong = cos( radpdg * rot ) |
---|
175 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
176 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
177 | else |
---|
178 | !* Normal case. Meteorological direction of wind related to true north. |
---|
179 | xpolg = strcmp(6) - strcmp(1) * xi0 |
---|
180 | ypolg = strcmp(5) - strcmp(1) * eta0 |
---|
181 | temp = sqrt ( xpolg ** 2 + ypolg ** 2 ) |
---|
182 | xpolg = xpolg / temp |
---|
183 | ypolg = ypolg / temp |
---|
184 | end if |
---|
185 | ue = ( ypolg * ug - xpolg * vg ) |
---|
186 | vn = ( ypolg * vg + xpolg * ug ) |
---|
187 | return |
---|
188 | end subroutine cg2cxy |
---|
189 | |
---|
190 | real function cgszll (strcmp, xlat,xlong) |
---|
191 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
192 | |
---|
193 | use par_mod, only: dp |
---|
194 | |
---|
195 | implicit none |
---|
196 | |
---|
197 | real :: strcmp(9), xlat, xlong |
---|
198 | |
---|
199 | real(kind=dp) :: slat,ymerc,efact |
---|
200 | |
---|
201 | if (xlat .gt. 89.985) then |
---|
202 | !* Close to north pole |
---|
203 | if (strcmp(1) .gt. 0.9999) then |
---|
204 | !* and to gamma == 1. |
---|
205 | cgszll = 2. * strcmp(7) |
---|
206 | return |
---|
207 | endif |
---|
208 | efact = cos(radpdg * xlat) |
---|
209 | if (efact .le. 0.) then |
---|
210 | cgszll = 0. |
---|
211 | return |
---|
212 | else |
---|
213 | ymerc = - log( efact /(1. + sin(radpdg * xlat))) |
---|
214 | endif |
---|
215 | else if (xlat .lt. -89.985) then |
---|
216 | !* Close to south pole |
---|
217 | if (strcmp(1) .lt. -0.9999) then |
---|
218 | !* and to gamma == -1.0 |
---|
219 | cgszll = 2. * strcmp(7) |
---|
220 | return |
---|
221 | endif |
---|
222 | efact = cos(radpdg * xlat) |
---|
223 | if (efact .le. 0.) then |
---|
224 | cgszll = 0. |
---|
225 | return |
---|
226 | else |
---|
227 | ymerc = log( efact /(1. - sin(radpdg * xlat))) |
---|
228 | endif |
---|
229 | else |
---|
230 | slat = sin(radpdg * xlat) |
---|
231 | ymerc = log((1. + slat) / (1. - slat))/2. |
---|
232 | !efact = exp(ymerc) |
---|
233 | !cgszll = 2. * strcmp(7) * exp (strcmp(1) * ymerc) |
---|
234 | !c / (efact + 1./efact) |
---|
235 | endif |
---|
236 | cgszll = strcmp(7) * cos(radpdg * xlat) * exp(strcmp(1) *ymerc) |
---|
237 | return |
---|
238 | end function cgszll |
---|
239 | |
---|
240 | real function cgszxy (strcmp, x,y) |
---|
241 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
242 | |
---|
243 | use par_mod, only: dp |
---|
244 | |
---|
245 | implicit none |
---|
246 | |
---|
247 | real :: strcmp(9) , x, y |
---|
248 | real(kind=dp) :: ymerc,efact, radial, temp |
---|
249 | real(kind=dp) :: xi0,eta0,xi,eta |
---|
250 | |
---|
251 | |
---|
252 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
253 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
254 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
255 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
256 | radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta) |
---|
257 | efact = strcmp(1) * radial |
---|
258 | if (efact .gt. almst1) then |
---|
259 | if (strcmp(1).gt.almst1) then |
---|
260 | cgszxy = 2. * strcmp(7) |
---|
261 | else |
---|
262 | cgszxy = 0. |
---|
263 | endif |
---|
264 | return |
---|
265 | endif |
---|
266 | if (abs(efact) .lt. 1.e-2) then |
---|
267 | temp = (efact / (2. - efact) )**2 |
---|
268 | ymerc = radial / (2. - efact) * (1. + temp * & |
---|
269 | (1./3. + temp * & |
---|
270 | (1./5. + temp * & |
---|
271 | (1./7. )))) |
---|
272 | else |
---|
273 | ymerc = - log( 1. - efact ) /2. /strcmp(1) |
---|
274 | endif |
---|
275 | if (ymerc .gt. 6.) then |
---|
276 | if (strcmp(1) .gt. almst1) then |
---|
277 | cgszxy = 2. * strcmp(7) |
---|
278 | else |
---|
279 | cgszxy = 0. |
---|
280 | endif |
---|
281 | else if (ymerc .lt. -6.) then |
---|
282 | if (strcmp(1) .lt. -almst1) then |
---|
283 | cgszxy = 2. * strcmp(7) |
---|
284 | else |
---|
285 | cgszxy = 0. |
---|
286 | endif |
---|
287 | else |
---|
288 | efact = exp(ymerc) |
---|
289 | cgszxy = 2. * strcmp(7) * exp (strcmp(1) * ymerc) & |
---|
290 | / (efact + 1./efact) |
---|
291 | endif |
---|
292 | return |
---|
293 | end function cgszxy |
---|
294 | |
---|
295 | subroutine cll2xy (strcmp, xlat,xlong, x,y) |
---|
296 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
297 | |
---|
298 | implicit none |
---|
299 | |
---|
300 | real :: strcmp(9) , xlat, xlong, x, y, xi, eta |
---|
301 | |
---|
302 | call cnllxy(strcmp, xlat,xlong, xi,eta) |
---|
303 | x = strcmp(3) + rearth/strcmp(7) * & |
---|
304 | (xi * strcmp(5) + eta * strcmp(6) ) |
---|
305 | y = strcmp(4) + rearth/strcmp(7) * & |
---|
306 | (eta * strcmp(5) - xi * strcmp(6) ) |
---|
307 | return |
---|
308 | end subroutine cll2xy |
---|
309 | |
---|
310 | subroutine cnllxy (strcmp, xlat,xlong, xi,eta) |
---|
311 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
312 | ! main transformation routine from latitude-longitude to |
---|
313 | ! canonical (equator-centered, radian unit) coordinates |
---|
314 | |
---|
315 | use par_mod, only: dp |
---|
316 | |
---|
317 | implicit none |
---|
318 | |
---|
319 | real :: strcmp(9), xlat, xlong, xi, eta, & |
---|
320 | gdlong, sndgam, csdgam, rhog1 |
---|
321 | real(kind=dp) :: gamma |
---|
322 | real(kind=dp) :: dlong,dlat,slat,mercy,gmercy |
---|
323 | |
---|
324 | gamma = strcmp(1) |
---|
325 | dlat = xlat |
---|
326 | dlong = cspanf(xlong - strcmp(2), -180., 180.) |
---|
327 | dlong = dlong * radpdg |
---|
328 | gdlong = gamma * dlong |
---|
329 | if (abs(gdlong) .lt. .01) then |
---|
330 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
331 | ! by zero in the case of mercator or near-mercator projections. |
---|
332 | gdlong = gdlong * gdlong |
---|
333 | sndgam = dlong * (1. - 1./6. * gdlong * & |
---|
334 | (1. - 1./20. * gdlong * & |
---|
335 | (1. - 1./42. * gdlong ))) |
---|
336 | csdgam = dlong * dlong * .5 * & |
---|
337 | (1. - 1./12. * gdlong * & |
---|
338 | (1. - 1./30. * gdlong * & |
---|
339 | (1. - 1./56. * gdlong ))) |
---|
340 | else |
---|
341 | ! Code for moderate values of gamma |
---|
342 | sndgam = sin (gdlong) /gamma |
---|
343 | csdgam = (1. - cos(gdlong) )/gamma /gamma |
---|
344 | endif |
---|
345 | slat = sin(radpdg * dlat) |
---|
346 | if ((slat .ge. almst1) .or. (slat .le. -almst1)) then |
---|
347 | eta = 1./strcmp(1) |
---|
348 | xi = 0. |
---|
349 | return |
---|
350 | endif |
---|
351 | mercy = .5 * log( (1. + slat) / (1. - slat) ) |
---|
352 | gmercy = gamma * mercy |
---|
353 | if (abs(gmercy) .lt. .001) then |
---|
354 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
355 | ! by zero in the case of mercator or near-mercator projections. |
---|
356 | rhog1 = mercy * (1. - .5 * gmercy * & |
---|
357 | (1. - 1./3. * gmercy * & |
---|
358 | (1. - 1./4. * gmercy ) ) ) |
---|
359 | else |
---|
360 | ! Code for moderate values of gamma |
---|
361 | rhog1 = (1. - exp(-gmercy)) / gamma |
---|
362 | endif |
---|
363 | eta = rhog1 + (1. - gamma * rhog1) * gamma * csdgam |
---|
364 | xi = (1. - gamma * rhog1 ) * sndgam |
---|
365 | end subroutine cnllxy |
---|
366 | |
---|
367 | subroutine cnxyll (strcmp, xi,eta, xlat,xlong) |
---|
368 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
369 | ! main transformation routine from canonical (equator-centered, |
---|
370 | ! radian unit) coordinates |
---|
371 | |
---|
372 | use par_mod, only: dp |
---|
373 | |
---|
374 | implicit none |
---|
375 | |
---|
376 | real :: strcmp(9), xlat, xlong, odist |
---|
377 | real(kind=dp) :: gamma,temp,arg1,arg2,ymerc,along,gxi,cgeta |
---|
378 | real(kind=dp) :: xi,eta |
---|
379 | |
---|
380 | gamma = strcmp(1) |
---|
381 | ! Calculate equivalent mercator coordinate |
---|
382 | odist = xi*xi + eta*eta |
---|
383 | arg2 = 2. * eta - gamma * (xi*xi + eta*eta) |
---|
384 | arg1 = gamma * arg2 |
---|
385 | ! Change by A. Stohl to avoid problems close to the poles |
---|
386 | ! if (arg1 .ge. almst1) then |
---|
387 | ! distance to north (or south) pole is zero (or imaginary ;) ) |
---|
388 | ! xlat = sign(90.,strcmp(1)) |
---|
389 | ! xlong = strcmp(2) |
---|
390 | ! return |
---|
391 | ! endif |
---|
392 | if (abs(arg1) .lt. .01) then |
---|
393 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
394 | ! by zero in the case of mercator or near-mercator projections. |
---|
395 | temp = (arg1 / (2. - arg1) )**2 |
---|
396 | ymerc = arg2 / (2. - arg1) * (1. + temp * & |
---|
397 | (1./3. + temp * & |
---|
398 | (1./5. + temp * & |
---|
399 | (1./7. )))) |
---|
400 | else |
---|
401 | ! Code for moderate values of gamma |
---|
402 | ymerc = - log ( 1. - arg1 ) /2. / gamma |
---|
403 | endif |
---|
404 | ! Convert ymerc to latitude |
---|
405 | temp = exp( - abs(ymerc) ) |
---|
406 | xlat = sign(atan2((1. - temp) * (1. + temp), 2. * temp), ymerc) |
---|
407 | ! Find longitudes |
---|
408 | gxi = gamma*xi |
---|
409 | cgeta = 1. - gamma * eta |
---|
410 | if ( abs(gxi) .lt. .01*cgeta ) then |
---|
411 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
412 | ! by zero in the case of mercator or near-mercator projections. |
---|
413 | temp = ( gxi /cgeta )**2 |
---|
414 | along = xi / cgeta * (1. - temp * & |
---|
415 | (1./3. - temp * & |
---|
416 | (1./5. - temp * & |
---|
417 | (1./7. )))) |
---|
418 | else |
---|
419 | ! Code for moderate values of gamma |
---|
420 | along = atan2( gxi, cgeta) / gamma |
---|
421 | endif |
---|
422 | xlong = sngl(strcmp(2) + dgprad * along) |
---|
423 | xlat = xlat * dgprad |
---|
424 | return |
---|
425 | end subroutine cnxyll |
---|
426 | |
---|
427 | subroutine cpolll (strcmp, xlat,xlong, enx,eny,enz) |
---|
428 | !* Written on 11/23/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
429 | |
---|
430 | use par_mod, only: dp |
---|
431 | |
---|
432 | implicit none |
---|
433 | |
---|
434 | real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot |
---|
435 | real :: strcmp(9), xlat, xlong, enx, eny, enz, clat |
---|
436 | |
---|
437 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
438 | rot = - strcmp(1) * along |
---|
439 | slong = sin( radpdg * rot ) |
---|
440 | clong = cos( radpdg * rot ) |
---|
441 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
442 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
443 | clat = cos(radpdg * xlat) |
---|
444 | enx = clat * xpolg |
---|
445 | eny = clat * ypolg |
---|
446 | enz = sin(radpdg * xlat) |
---|
447 | return |
---|
448 | end subroutine cpolll |
---|
449 | |
---|
450 | subroutine cpolxy (strcmp, x,y, enx,eny,enz) |
---|
451 | !* Written on 11/26/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
452 | |
---|
453 | use par_mod, only: dp |
---|
454 | |
---|
455 | implicit none |
---|
456 | |
---|
457 | real :: strcmp(9) , x, y, enx, eny, enz |
---|
458 | real(kind=dp) :: xpol,ypol,temp,xi0,eta0,xi,eta,radial |
---|
459 | real(kind=dp) :: temp2,ymerc,arg,oarg,clat |
---|
460 | |
---|
461 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
462 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
463 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
464 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
465 | radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta) |
---|
466 | temp = strcmp(1) * radial |
---|
467 | if (temp .ge. 1.) then |
---|
468 | enx = 0. |
---|
469 | eny = 0. |
---|
470 | enz = sign(1.,strcmp(1)) |
---|
471 | return |
---|
472 | endif |
---|
473 | if (abs(temp).lt.1.e-2) then |
---|
474 | temp2 = (temp / (2. - temp))**2 |
---|
475 | ymerc = radial / (2. - temp) * (1. + temp2 * & |
---|
476 | (1./3. + temp2 * & |
---|
477 | (1./5. + temp2 * & |
---|
478 | (1./7.)))) |
---|
479 | else |
---|
480 | ymerc = -.5 * log(1. - temp) / strcmp(1) |
---|
481 | endif |
---|
482 | arg = exp( ymerc ) |
---|
483 | oarg = 1./arg |
---|
484 | clat = 2./(arg + oarg) |
---|
485 | enz = (arg - oarg) * clat /2. |
---|
486 | temp = clat / sqrt(1. - temp) |
---|
487 | xpol = - xi * strcmp(1) * temp |
---|
488 | ypol = (1. - eta * strcmp(1) ) * temp |
---|
489 | enx = xpol * strcmp(5) + ypol * strcmp(6) |
---|
490 | eny = ypol * strcmp(5) - xpol * strcmp(6) |
---|
491 | return |
---|
492 | end subroutine cpolxy |
---|
493 | |
---|
494 | real function cspanf (value, begin, end) |
---|
495 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
496 | !* real function cspanf returns a value in the interval (begin,end] |
---|
497 | !* which is equivalent to value, mod (end - begin). It is used to |
---|
498 | !* reduce periodic variables to a standard range. It adjusts for the |
---|
499 | !* behavior of the mod function which provides positive results for |
---|
500 | !* positive input, and negative results for negative input |
---|
501 | !* input: |
---|
502 | !* value - real number to be reduced to the span |
---|
503 | !* begin - first value of the span |
---|
504 | !* end - last value of the span |
---|
505 | !* returns: |
---|
506 | !* the reduced value |
---|
507 | !* examples: |
---|
508 | !* along = cspanf(xlong, -180., +180.) |
---|
509 | !* dir = cspanf(angle, 0., 360.) |
---|
510 | |
---|
511 | implicit none |
---|
512 | |
---|
513 | real :: first,last, value, begin, end, val |
---|
514 | |
---|
515 | first = min(begin,end) |
---|
516 | last = max(begin,end) |
---|
517 | val = mod( value - first , last - first) |
---|
518 | if ( val .le. 0.) then |
---|
519 | cspanf = val + last |
---|
520 | else |
---|
521 | cspanf = val + first |
---|
522 | endif |
---|
523 | return |
---|
524 | end function cspanf |
---|
525 | |
---|
526 | subroutine cxy2ll (strcmp, x,y, xlat,xlong) |
---|
527 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
528 | |
---|
529 | use par_mod, only: dp |
---|
530 | |
---|
531 | implicit none |
---|
532 | |
---|
533 | real(kind=dp) :: xi0,eta0,xi,eta |
---|
534 | real :: strcmp(9), x, y, xlat, xlong |
---|
535 | |
---|
536 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
537 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
538 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
539 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
540 | call cnxyll(strcmp, xi,eta, xlat,xlong) |
---|
541 | xlong = cspanf(xlong, -180., 180.) |
---|
542 | return |
---|
543 | end subroutine cxy2ll |
---|
544 | |
---|
545 | real function eqvlat (xlat1,xlat2) |
---|
546 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
547 | |
---|
548 | implicit none |
---|
549 | |
---|
550 | real :: xlat1, xlat2, x, ssind, sinl1, sinl2, al1, al2, tau |
---|
551 | |
---|
552 | ssind(x) = sin (radpdg*x) |
---|
553 | sinl1 = ssind (xlat1) |
---|
554 | sinl2 = ssind (xlat2) |
---|
555 | if (abs(sinl1 - sinl2) .gt. .001) then |
---|
556 | al1 = log((1. - sinl1)/(1. - sinl2)) |
---|
557 | al2 = log((1. + sinl1)/(1. + sinl2)) |
---|
558 | else |
---|
559 | ! Case lat1 near or equal to lat2 |
---|
560 | tau = - (sinl1 - sinl2)/(2. - sinl1 - sinl2) |
---|
561 | tau = tau*tau |
---|
562 | al1 = 2. / (2. - sinl1 - sinl2) * (1. + tau * & |
---|
563 | (1./3. + tau * & |
---|
564 | (1./5. + tau * & |
---|
565 | (1./7.)))) |
---|
566 | tau = (sinl1 - sinl2)/(2. + sinl1 + sinl2) |
---|
567 | tau = tau*tau |
---|
568 | al2 = -2. / (2. + sinl1 + sinl2) * (1. + tau * & |
---|
569 | (1./3. + tau * & |
---|
570 | (1./5. + tau * & |
---|
571 | (1./7.)))) |
---|
572 | endif |
---|
573 | eqvlat = asin((al1 + al2) / (al1 - al2))/radpdg |
---|
574 | return |
---|
575 | end function eqvlat |
---|
576 | |
---|
577 | subroutine stcm1p(strcmp, x1,y1, xlat1,xlong1, & |
---|
578 | xlatg,xlongg, gridsz, orient) |
---|
579 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
580 | |
---|
581 | implicit none |
---|
582 | |
---|
583 | integer :: k |
---|
584 | real :: strcmp(9), x1, y1, xlat1, xlong1, turn, orient, & |
---|
585 | xlatg, xlongg, gridsz, x1a, y1a |
---|
586 | |
---|
587 | do k=3,4 |
---|
588 | strcmp (k) = 0. |
---|
589 | enddo |
---|
590 | turn = radpdg * (orient - strcmp(1) * & |
---|
591 | cspanf(xlongg - strcmp(2), -180., 180.) ) |
---|
592 | strcmp (5) = cos (turn) |
---|
593 | strcmp (6) = - sin (turn) |
---|
594 | strcmp (7) = 1. |
---|
595 | strcmp (7) = gridsz * strcmp(7) & |
---|
596 | / cgszll(strcmp, xlatg, strcmp(2)) |
---|
597 | call cll2xy (strcmp, xlat1,xlong1, x1a,y1a) |
---|
598 | strcmp(3) = strcmp(3) + x1 - x1a |
---|
599 | strcmp(4) = strcmp(4) + y1 - y1a |
---|
600 | return |
---|
601 | end subroutine stcm1p |
---|
602 | |
---|
603 | subroutine stcm2p(strcmp, x1,y1, xlat1,xlong1, & |
---|
604 | x2,y2, xlat2,xlong2) |
---|
605 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
606 | |
---|
607 | implicit none |
---|
608 | |
---|
609 | real :: strcmp(9), x1, y1, xlat1, xlong1, & |
---|
610 | x2, y2, xlat2, xlong2 |
---|
611 | |
---|
612 | integer :: k |
---|
613 | real :: x1a, y1a, x2a, y2a, den, dena |
---|
614 | |
---|
615 | do k=3,6 |
---|
616 | strcmp (k) = 0. |
---|
617 | enddo |
---|
618 | strcmp (5) = 1. |
---|
619 | strcmp (7) = 1. |
---|
620 | call cll2xy (strcmp, xlat1,xlong1, x1a,y1a) |
---|
621 | call cll2xy (strcmp, xlat2,xlong2, x2a,y2a) |
---|
622 | den = sqrt( (x1 - x2)**2 + (y1 - y2)**2 ) |
---|
623 | dena = sqrt( (x1a - x2a)**2 + (y1a - y2a)**2 ) |
---|
624 | strcmp(5) = ((x1a - x2a)*(x1 - x2) + (y1a - y2a) * (y1 - y2)) & |
---|
625 | /den /dena |
---|
626 | strcmp(6) = ((y1a - y2a)*(x1 - x2) - (x1a - x2a) * (y1 - y2)) & |
---|
627 | /den /dena |
---|
628 | strcmp (7) = strcmp(7) * dena / den |
---|
629 | call cll2xy (strcmp, xlat1,xlong1, x1a,y1a) |
---|
630 | strcmp(3) = strcmp(3) + x1 - x1a |
---|
631 | strcmp(4) = strcmp(4) + y1 - y1a |
---|
632 | return |
---|
633 | end subroutine stcm2p |
---|
634 | |
---|
635 | !* General conformal map routines for meteorological modelers |
---|
636 | !* written on 3/31/94 by |
---|
637 | |
---|
638 | !* Dr. Albion Taylor |
---|
639 | !* NOAA / OAR / ARL phone: (301) 713-0295 x 132 |
---|
640 | !* rm. 3151, 1315 east-west highway fax: (301) 713-0119 |
---|
641 | !* silver spring, md 20910 e-mail: adtaylor@arlrisc.ssmc.noaa.gov |
---|
642 | |
---|
643 | !* subroutine stlmbr (strcmp, tnglat, clong) |
---|
644 | !* This routine initializes the map structure array strcmp to |
---|
645 | !* the form of a specific map projection |
---|
646 | !* inputs: |
---|
647 | !* tnglat - the latitude at which the projection will be tangent |
---|
648 | !* to the earth. +90. For north polar stereographic, |
---|
649 | !* -90. for south polar stereographic, 0. For mercator, |
---|
650 | !* and other values for lambert conformal. |
---|
651 | !* -90 <= tnglat <= 90. |
---|
652 | !* clong - a longitude in the region under consideration. Longitudes |
---|
653 | !* between clong-180. and clong+180. Will be mapped in one |
---|
654 | !* connected region |
---|
655 | !* outputs: |
---|
656 | !* strcmp - a 9-value map structure array for use with subsequent |
---|
657 | !* calls to the coordinate transform routines. |
---|
658 | !* |
---|
659 | !* real function eqvlat (xlat1,xlat2) |
---|
660 | !* This function is provided to assist in finding the tangent latitude |
---|
661 | !* equivalent to the 2-reference latitude specification in the legend |
---|
662 | !* of most lambert conformal maps. If the map specifies "scale |
---|
663 | !* 1:xxxxx true at 40n and 60n", then eqvlat(40.,60.) will return the |
---|
664 | !* equivalent tangent latitude. |
---|
665 | !* inputs: |
---|
666 | !* xlat1,xlat2: the two latitudes specified in the map legend |
---|
667 | !* returns: |
---|
668 | !* the equivalent tangent latitude |
---|
669 | !* example: call stlmbr(strcmp, eqvlat(40.,60.), 90.) |
---|
670 | |
---|
671 | !* subroutine stcm2p (strcmp, x1,y1, xlat1,xlong1, |
---|
672 | !* x2,y2, xlat2,xlong2) |
---|
673 | !* subroutine stcm1p (strcmp, x1,y1, xlat1,xlong1, |
---|
674 | !* xlatg,xlongg, gridsz, orient) |
---|
675 | !* These routines complete the specification of the map structure |
---|
676 | !* array by conforming the map coordinates to the specifications |
---|
677 | !* of a particular grid. Either stcm1p or stcm2p must be called, |
---|
678 | !* but not both |
---|
679 | !* inputs: |
---|
680 | !* strcmp - a 9-value map structure array, set to a particular map |
---|
681 | !* form by a previous call to stlmbr |
---|
682 | !* for stcm2p: |
---|
683 | !* x1,y1, x2,y2 - the map coordinates of two points on the grid |
---|
684 | !* xlat1,xlong1, xlat2,xlong2 - the geographic coordinates of the |
---|
685 | !* same two points |
---|
686 | !* for stcm1p: |
---|
687 | !* x1,y1 - the map coordinates of one point on the grid |
---|
688 | !* xlat1,xlong1 - the geographic coordinates of the same point |
---|
689 | !* xlatg,xlongg - latitude and longitude of reference point for |
---|
690 | !* gridsz and orientation specification. |
---|
691 | !* gridsz - the desired grid size in kilometers, at xlatg,xlongg |
---|
692 | !* orient - the angle, with respect to north, of a y-grid line, at |
---|
693 | !* the point xlatg,xlongg |
---|
694 | !* outputs: |
---|
695 | !* strcmp - a 9-value map structure array, fully set for use by |
---|
696 | !* other subroutines in this system |
---|
697 | |
---|
698 | !* subroutine cll2xy (strcmp, xlat,xlong, x,y) |
---|
699 | !* subroutine cxy2ll (strcmp, x,y, xlat,xlong) |
---|
700 | !* these routines convert between map coordinates x,y |
---|
701 | !* and geographic coordinates xlat,xlong |
---|
702 | !* inputs: |
---|
703 | !* strcmp(9) - 9-value map structure array |
---|
704 | !* for cll2xy: xlat,xlong - geographic coordinates |
---|
705 | !* for cxy2ll: x,y - map coordinates |
---|
706 | !* outputs: |
---|
707 | !* for cll2xy: x,y - map coordinates |
---|
708 | !* for cxy2ll: xlat,xlong - geographic coordinates |
---|
709 | |
---|
710 | !* subroutine cc2gxy (strcmp, x,y, ue,vn, ug,vg) |
---|
711 | !* subroutine cg2cxy (strcmp, x,y, ug,vg, ue,vn) |
---|
712 | !* subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg) |
---|
713 | !* subroutine cg2cll (strcmp, xlat,xlong, ug,vg, ue,vn) |
---|
714 | !* These subroutines convert vector wind components from |
---|
715 | !* geographic, or compass, coordinates, to map or |
---|
716 | !* grid coordinates. The site of the wind to be |
---|
717 | !* converted may be given either in geographic or |
---|
718 | !* map coordinates. Wind components are all in kilometers |
---|
719 | !* per hour, whether geographic or map coordinates. |
---|
720 | !* inputs: |
---|
721 | !* strcmp(9) - 9-value map structure array |
---|
722 | !* for cc2gxy and cg2cxy: x,y - map coordinates of site |
---|
723 | !* for cc2gll and cg2cll: xlat,xlong - geographic coordinates of site |
---|
724 | !* for cc2gxy and cc2gll: ue,vn - east and north wind components |
---|
725 | !* for cg2cxy and cg2cll: ug,vg - x- and y- direction wind components |
---|
726 | !* outputs: |
---|
727 | !* for cc2gxy and cc2gll: ug,vg - x- and y- direction wind components |
---|
728 | !* for cg2cxy and cg2cll: ue,vn - east and north wind components |
---|
729 | |
---|
730 | !* subroutine ccrvxy (strcmp, x, y, gx,gy) |
---|
731 | !* subroutine ccrvll (strcmp, xlat,xlong, gx,gy) |
---|
732 | !* These subroutines return the curvature vector (gx,gy), as referenced |
---|
733 | !* to map coordinates, induced by the map transformation. When |
---|
734 | !* non-linear terms in wind speed are important, a "geodesic" force |
---|
735 | !* should be included in the vector form [ (u,u) g - (u,g) u ] where the |
---|
736 | !* inner product (u,g) is defined as ux*gx + uy*gy. |
---|
737 | !* inputs: |
---|
738 | !* strcmp(9) - 9-value map structure array |
---|
739 | !* for ccrvxy: x,y - map coordinates of site |
---|
740 | !* for ccrvll: xlat,xlong - geographic coordinates of site |
---|
741 | !* outputs: |
---|
742 | !* gx,gy - vector coefficients of curvature, in units radians |
---|
743 | !* per kilometer |
---|
744 | |
---|
745 | !* real function cgszll (strcmp, xlat,xlong) |
---|
746 | !* real function cgszxy (strcmp, x,y) |
---|
747 | !* These functions return the size, in kilometers, of each unit of |
---|
748 | !* motion in map coordinates (grid size). The grid size at any |
---|
749 | !* location depends on that location; the position may be given in |
---|
750 | !* either map or geographic coordinates. |
---|
751 | !* inputs: |
---|
752 | !* strcmp(9) - 9-value map structure array |
---|
753 | !* for cgszxy: x,y - map coordinates of site |
---|
754 | !* for cgszll: xlat,xlong - geographic coordinates of site |
---|
755 | !* returns: |
---|
756 | !* gridsize in kilometers at given site. |
---|
757 | |
---|
758 | !* subroutine cpolxy (strcmp, x,y, enx,eny,enz) |
---|
759 | !* subroutine cpolll (strcmp, xlat,xlong, enx,eny,enz) |
---|
760 | !* These subroutines provide 3-d vector components of a unit vector |
---|
761 | !* in the direction of the north polar axis. When multiplied |
---|
762 | !* by twice the rotation rate of the earth (2 * pi/24 hr), the |
---|
763 | !* vertical component yields the coriolis factor. |
---|
764 | !* inputs: |
---|
765 | !* strcmp(9) - 9-value map structure array |
---|
766 | !* for cpolxy: x,y - map coordinates of site |
---|
767 | !* for cpolll: xlat,xlong - geographic coordinates of site |
---|
768 | !* returns: |
---|
769 | !* enx,eny,enz the direction cosines of a unit vector in the |
---|
770 | !* direction of the rotation axis of the earth |
---|
771 | |
---|
772 | !* subroutine cnllxy (strcmp, xlat,xlong, xi,eta) |
---|
773 | !* subroutine cnxyll (strcmp, xi,eta, xlat,xlong) |
---|
774 | !* These subroutines perform the underlying transformations from |
---|
775 | !* geographic coordinates to and from canonical (equator centered) |
---|
776 | !* coordinates. They are called by cxy2ll and cll2xy, but are not |
---|
777 | !* intended to be called directly |
---|
778 | |
---|
779 | !* real function cspanf (value, begin, end) |
---|
780 | !* This function assists other routines in providing a longitude in |
---|
781 | !* the proper range. It adds to value whatever multiple of |
---|
782 | !* (end - begin) is needed to return a number begin < cspanf <= end |
---|
783 | |
---|
784 | subroutine stlmbr(strcmp, tnglat, xlong) |
---|
785 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
786 | |
---|
787 | implicit none |
---|
788 | |
---|
789 | real :: strcmp(9), tnglat, xlong |
---|
790 | |
---|
791 | real :: eta, xi |
---|
792 | |
---|
793 | strcmp(1) = sin(radpdg * tnglat) |
---|
794 | !* gamma = sine of the tangent latitude |
---|
795 | strcmp(2) = cspanf( xlong, -180., +180.) |
---|
796 | !* lambda_0 = reference longitude |
---|
797 | strcmp(3) = 0. |
---|
798 | !* x_0 = x- grid coordinate of origin (xi,eta) = (0.,0.) |
---|
799 | strcmp(4) = 0. |
---|
800 | !* y_0 = y-grid coordinate of origin (xi,eta) = (0.,0.) |
---|
801 | strcmp(5) = 1. |
---|
802 | !* Cosine of rotation angle from xi,eta to x,y |
---|
803 | strcmp(6) = 0. |
---|
804 | !* Sine of rotation angle from xi,eta to x,y |
---|
805 | strcmp(7) = rearth |
---|
806 | !* Gridsize in kilometers at the equator |
---|
807 | call cnllxy(strcmp, 89.,xlong, xi,eta) |
---|
808 | strcmp(8) = 2. * eta - strcmp(1) * eta * eta |
---|
809 | !* Radial coordinate for 1 degree from north pole |
---|
810 | call cnllxy(strcmp, -89.,xlong, xi,eta) |
---|
811 | strcmp(9) = 2. * eta - strcmp(1) * eta * eta |
---|
812 | !* Radial coordinate for 1 degree from south pole |
---|
813 | return |
---|
814 | end subroutine stlmbr |
---|
815 | |
---|
816 | end module cmapf_mod |
---|