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