35 real,
parameter :: rearth=6371.2, almst1=.9999999
37 real,
parameter :: pi=3.14159265358979
38 real,
parameter :: radpdg=pi/180., dgprad=180./pi
42 subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg)
49 real :: strcmp(9), xlat, xlong, ue, vn, ug, vg
51 real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot
53 along =
cspanf( xlong - strcmp(2), -180., 180.)
54 if (xlat.gt.89.985)
then
56 rot = - strcmp(1) * along + xlong - 180.
57 elseif (xlat.lt.-89.985)
then
59 rot = - strcmp(1) * along - xlong
61 rot = - strcmp(1) * along
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
72 subroutine ccrvll (strcmp, xlat,xlong, gx,gy)
79 real(kind=dp) :: xpolg,ypolg,temp,along,slong,clong,ctemp, curv
80 real :: strcmp(9), xlat, xlong, gx, gy
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
95 subroutine ccrvxy (strcmp, x,y, gx,gy)
102 real :: strcmp(9), x, y, gx, gy
103 real(kind=dp) :: xpolg,ypolg,temp,ymerc,efact,curv
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 )
110 ymerc = - log( temp) /strcmp(1)
112 curv = ( (strcmp(1) - 1.d0) * efact + &
113 (strcmp(1) + 1.d0) / efact ) &
115 gx = xpolg * curv / temp
116 gy = ypolg * curv / temp
118 if (abs(strcmp(1)) .eq. 1.)
then
129 subroutine cg2cll (strcmp, xlat,xlong, ug,vg, ue,vn)
136 real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot
137 real :: strcmp(9), xlat, xlong, ug, vg, ue, vn
139 along =
cspanf( xlong - strcmp(2), -180., 180.)
140 if (xlat.gt.89.985)
then
142 rot = - strcmp(1) * along + xlong - 180.
143 elseif (xlat.lt.-89.985)
then
145 rot = - strcmp(1) * along - xlong
147 rot = - strcmp(1) * along
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
158 subroutine cg2cxy (strcmp, x,y, ug,vg, ue,vn)
165 real :: strcmp(9) , x, y, ug, vg, ue, vn
167 real :: clong, radial, rot, slong, xlat, xlong
168 real(kind=dp) :: xpolg,ypolg,temp,xi0,eta0,xi,eta
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
178 call
cnxyll(strcmp, xi,eta, xlat,xlong)
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
188 call
cnxyll(strcmp, xi,eta, xlat,xlong)
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)
197 xpolg = strcmp(6) - strcmp(1) * xi0
198 ypolg = strcmp(5) - strcmp(1) * eta0
199 temp = sqrt( xpolg ** 2 + ypolg ** 2 )
203 ue = ( ypolg * ug - xpolg * vg )
204 vn = ( ypolg * vg + xpolg * ug )
208 real function cgszll (strcmp, xlat,xlong)
215 real :: strcmp(9), xlat, xlong
217 real(kind=dp) :: slat,ymerc,efact
219 if (xlat .gt. 89.985)
then
221 if (strcmp(1) .gt. 0.9999)
then
226 efact = cos(radpdg * xlat)
227 if (efact .le. 0.)
then
231 ymerc = - log( efact /(1. + sin(radpdg * xlat)))
233 else if (xlat .lt. -89.985)
then
235 if (strcmp(1) .lt. -0.9999)
then
240 efact = cos(radpdg * xlat)
241 if (efact .le. 0.)
then
245 ymerc = log( efact /(1. - sin(radpdg * xlat)))
248 slat = sin(radpdg * xlat)
249 ymerc = log((1. + slat) / (1. - slat))/2.
254 cgszll = strcmp(7) * cos(radpdg * xlat) * exp(strcmp(1) *ymerc)
265 real :: strcmp(9) , x, y
266 real(kind=dp) :: ymerc,efact, radial, temp
267 real(kind=dp) :: xi0,eta0,xi,eta
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
284 if (abs(efact) .lt. 1.e-2)
then
285 temp = (efact / (2. - efact) )**2
286 ymerc = radial / (2. - efact) * (1. + temp * &
291 ymerc = - log( 1. - efact ) /2. /strcmp(1)
293 if (ymerc .gt. 6.)
then
294 if (strcmp(1) .gt. almst1)
then
299 else if (ymerc .lt. -6.)
then
300 if (strcmp(1) .lt. -almst1)
then
307 cgszxy = 2. * strcmp(7) * exp(strcmp(1) * ymerc) &
313 subroutine cll2xy (strcmp, xlat,xlong, x,y)
318 real :: strcmp(9) , xlat, xlong, x, y, xi, eta
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) )
328 subroutine cnllxy (strcmp, xlat,xlong, xi,eta)
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
344 dlong =
cspanf(xlong - strcmp(2), -180., 180.)
345 dlong = dlong * radpdg
346 gdlong = gamma * dlong
347 if (abs(gdlong) .lt. .01)
then
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 )))
360 sndgam = sin(gdlong) /gamma
361 csdgam = (1. - cos(gdlong) )/gamma /gamma
363 slat = sin(radpdg * dlat)
364 if ((slat .ge. almst1) .or. (slat .le. -almst1))
then
369 mercy = .5 * log( (1. + slat) / (1. - slat) )
370 gmercy = gamma * mercy
371 if (abs(gmercy) .lt. .001)
then
374 rhog1 = mercy * (1. - .5 * gmercy * &
375 (1. - 1./3. * gmercy * &
376 (1. - 1./4. * gmercy ) ) )
379 rhog1 = (1. - exp(-gmercy)) / gamma
381 eta = rhog1 + (1. - gamma * rhog1) * gamma * csdgam
382 xi = (1. - gamma * rhog1 ) * sndgam
385 subroutine cnxyll (strcmp, xi,eta, xlat,xlong)
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
400 odist = xi*xi + eta*eta
401 arg2 = 2. * eta - gamma * (xi*xi + eta*eta)
410 if (abs(arg1) .lt. .01)
then
413 temp = (arg1 / (2. - arg1) )**2
414 ymerc = arg2 / (2. - arg1) * (1. + temp * &
420 ymerc = - log( 1. - arg1 ) /2. / gamma
423 temp = exp( - abs(ymerc) )
424 xlat = sign(atan2((1. - temp) * (1. + temp), 2. * temp), ymerc)
427 cgeta = 1. - gamma * eta
428 if ( abs(gxi) .lt. .01*cgeta )
then
431 temp = ( gxi /cgeta )**2
432 along = xi / cgeta * (1. - temp * &
438 along = atan2( gxi, cgeta) / gamma
440 xlong = sngl(strcmp(2) + dgprad * along)
445 subroutine cpolll (strcmp, xlat,xlong, enx,eny,enz)
452 real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot
453 real :: strcmp(9), xlat, xlong, enx, eny, enz, clat
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)
464 enz = sin(radpdg * xlat)
468 subroutine cpolxy (strcmp, x,y, enx,eny,enz)
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
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
488 enz = sign(1.,strcmp(1))
491 if (abs(temp).lt.1.e-2)
then
492 temp2 = (temp / (2. - temp))**2
493 ymerc = radial / (2. - temp) * (1. + temp2 * &
498 ymerc = -.5 * log(1. - temp) / strcmp(1)
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)
531 real :: first,last, value, begin, end, val
533 first = min(begin,end)
534 last = max(begin,end)
535 val = mod( value - first , last - first)
536 if ( val .le. 0.)
then
544 subroutine cxy2ll (strcmp, x,y, xlat,xlong)
551 real(kind=dp) :: xi0,eta0,xi,eta
552 real :: strcmp(9), x, y, xlat, xlong
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.)
568 real :: xlat1, xlat2, x, ssind, sinl1, sinl2, al1, al2, tau
570 ssind(x) = sin(radpdg*x)
573 if (abs(sinl1 - sinl2) .gt. .001)
then
574 al1 = log((1. - sinl1)/(1. - sinl2))
575 al2 = log((1. + sinl1)/(1. + sinl2))
578 tau = - (sinl1 - sinl2)/(2. - sinl1 - sinl2)
580 al1 = 2. / (2. - sinl1 - sinl2) * (1. + tau * &
584 tau = (sinl1 - sinl2)/(2. + sinl1 + sinl2)
586 al2 = -2. / (2. + sinl1 + sinl2) * (1. + tau * &
591 eqvlat = asin((al1 + al2) / (al1 - al2))/radpdg
595 subroutine stcm1p(strcmp, x1,y1, xlat1,xlong1, &
596 xlatg,xlongg, gridsz, orient)
602 real :: strcmp(9), x1, y1, xlat1, xlong1, turn, orient, &
603 xlatg, xlongg, gridsz, x1a, y1a
608 turn = radpdg * (orient - strcmp(1) * &
609 cspanf(xlongg - strcmp(2), -180., 180.) )
610 strcmp(5) = cos(turn)
611 strcmp(6) = - sin(turn)
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
621 subroutine stcm2p(strcmp, x1,y1, xlat1,xlong1, &
627 real :: strcmp(9), x1, y1, xlat1, xlong1, &
628 x2, y2, xlat2, xlong2
631 real :: x1a, y1a, x2a, y2a, den, dena
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)) &
644 strcmp(6) = ((y1a - y2a)*(x1 - x2) - (x1a - x2a) * (y1 - y2)) &
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
807 real :: strcmp(9), tnglat, xlong
811 strcmp(1) = sin(radpdg * tnglat)
813 strcmp(2) =
cspanf( xlong, -180., +180.)
825 call
cnllxy(strcmp, 89.,xlong, xi,eta)
826 strcmp(8) = 2. * eta - strcmp(1) * eta * eta
828 call
cnllxy(strcmp, -89.,xlong, xi,eta)
829 strcmp(9) = 2. * eta - strcmp(1) * eta * eta
real function, public cgszll(strcmp, xlat, xlong)
real function cspanf(value, begin, end)
subroutine stcm1p(strcmp, x1, y1, xlat1, xlong1, xlatg, xlongg, gridsz, orient)
subroutine, public stcm2p(strcmp, x1, y1, xlat1, xlong1, x2, y2, xlat2, xlong2)
subroutine, public cc2gll(strcmp, xlat, xlong, ue, vn, ug, vg)
subroutine cnxyll(strcmp, xi, eta, xlat, xlong)
subroutine, public cxy2ll(strcmp, x, y, xlat, xlong)
subroutine cg2cxy(strcmp, x, y, ug, vg, ue, vn)
subroutine cnllxy(strcmp, xlat, xlong, xi, eta)
real function eqvlat(xlat1, xlat2)
subroutine ccrvll(strcmp, xlat, xlong, gx, gy)
subroutine cpolll(strcmp, xlat, xlong, enx, eny, enz)
subroutine cpolxy(strcmp, x, y, enx, eny, enz)
subroutine, public stlmbr(strcmp, tnglat, xlong)
subroutine ccrvxy(strcmp, x, y, gx, gy)
real function cgszxy(strcmp, x, y)
subroutine, public cll2xy(strcmp, xlat, xlong, x, y)
subroutine cg2cll(strcmp, xlat, xlong, ug, vg, ue, vn)