source: flexpart.git/src/cmapf_mod.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 23547f3, checked in by Ignacio Pisso <ip@…>, 4 years ago

remove tabs from files of the form src/*.f90

  • Property mode set to 100644
File size: 25.7 KB
Line 
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
8module 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
22contains
23
24subroutine 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
52end subroutine cc2gll
53
54subroutine 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
75end subroutine ccrvll
76
77subroutine 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
109end subroutine ccrvxy
110
111subroutine 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
138end subroutine cg2cll
139
140subroutine 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
188end subroutine cg2cxy
189
190real 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
238end function cgszll
239
240real 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
293end function cgszxy
294
295subroutine 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
308end subroutine cll2xy
309
310subroutine 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
365end subroutine cnllxy
366
367subroutine 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
425end subroutine cnxyll
426
427subroutine 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
448end subroutine cpolll
449
450subroutine 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
492end subroutine cpolxy
493
494real 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
524end function cspanf
525
526subroutine 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
543end subroutine cxy2ll
544
545real 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
575end function eqvlat
576
577subroutine 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
601end subroutine stcm1p
602
603subroutine 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
633end 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
784subroutine 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
814end subroutine stlmbr
815
816end module cmapf_mod
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG