source: trunk/src/cmapf_mod.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 27.0 KB
Line 
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
26module 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
40contains
41
42subroutine 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
70end subroutine cc2gll
71
72subroutine 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
93end subroutine ccrvll
94
95subroutine 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
127end subroutine ccrvxy
128
129subroutine 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
156end subroutine cg2cll
157
158subroutine 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
206end subroutine cg2cxy
207
208real 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
256end function cgszll
257
258real 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
311end function cgszxy
312
313subroutine 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
326end subroutine cll2xy
327
328subroutine 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
383end subroutine cnllxy
384
385subroutine 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
443end subroutine cnxyll
444
445subroutine 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
466end subroutine cpolll
467
468subroutine 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
510end subroutine cpolxy
511
512real 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
542end function cspanf
543
544subroutine 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
561end subroutine cxy2ll
562
563real 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
593end function eqvlat
594
595subroutine 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
619end subroutine stcm1p
620
621subroutine 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
651end 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
802subroutine 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
832end subroutine stlmbr
833
834end module cmapf_mod
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG