cmapf_mod.f90 25.6 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
! Changes to the routines by A. Stohl
! xi,xi0,eta,eta0 are double precision variables to avoid problems
! at poles

module cmapf_mod

  use par_mod, only: dp

  implicit none
  private

  public :: cc2gll, cll2xy, cgszll, cxy2ll, stlmbr, stcm2p

  real,parameter :: rearth=6371.2, almst1=.9999999

  real,parameter :: pi=3.14159265358979
  real,parameter :: radpdg=pi/180., dgprad=180./pi

contains

subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real :: strcmp(9), xlat, xlong, ue, vn, ug, vg

  real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot

  along = cspanf( xlong - strcmp(2), -180., 180.)
  if (xlat.gt.89.985) then
  !*  North polar meteorological orientation: "north" along prime meridian
    rot = - strcmp(1) * along + xlong - 180.
  elseif (xlat.lt.-89.985) then
  !*  South polar meteorological orientation: "north" along prime meridian
    rot = - strcmp(1) * along - xlong
  else
    rot = - strcmp(1) * along
  endif
  slong = sin( radpdg * rot )
  clong = cos( radpdg * rot )
  xpolg = slong * strcmp(5) + clong * strcmp(6)
  ypolg = clong * strcmp(5) - slong * strcmp(6)
  ug = ypolg * ue + xpolg * vn
  vg = ypolg * vn - xpolg * ue
  return
end subroutine cc2gll

subroutine ccrvll (strcmp, xlat,xlong, gx,gy)
  !*  Written on 9/20/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real(kind=dp) :: xpolg,ypolg,temp,along,slong,clong,ctemp, curv
  real :: strcmp(9), xlat, xlong, gx, gy

  along = cspanf( xlong - strcmp(2), -180., 180.)
  slong = sin( radpdg * strcmp(1) * along)
  clong = cos( radpdg * strcmp(1) * along)
  xpolg = - slong * strcmp(5) + clong * strcmp(6)
  ypolg = clong * strcmp(5) + slong * strcmp(6)
  temp = sin(radpdg * xlat)
  ctemp = cos(radpdg * xlat)
  curv = (strcmp(1) - temp) / ctemp / rearth
  gx = curv * xpolg
  gy = curv * ypolg
  return
end subroutine ccrvll

subroutine ccrvxy (strcmp, x,y, gx,gy)
  !*  Written on 9/20/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real :: strcmp(9), x, y, gx, gy
  real(kind=dp) :: xpolg,ypolg,temp,ymerc,efact,curv

  temp = strcmp(1) * strcmp(7) /rearth
  xpolg = strcmp(6) + temp * (strcmp(3) - x)
  ypolg = strcmp(5) + temp * (strcmp(4) - y)
  temp = sqrt ( xpolg ** 2 + ypolg ** 2 )
  if (temp.gt.0.) then
    ymerc = - log( temp) /strcmp(1)
    efact = exp(ymerc)
    curv = ( (strcmp(1) - 1.d0) * efact + &
         (strcmp(1) + 1.d0) / efact ) &
         * .5d0 / rearth
    gx = xpolg * curv / temp
    gy = ypolg * curv / temp
  else
    if (abs(strcmp(1)) .eq. 1.) then
      gx = 0.
      gy = 0.
    else
      gx = 1./rearth
      gy = 1./rearth
    endif
  endif
  return
end subroutine ccrvxy

subroutine cg2cll (strcmp, xlat,xlong, ug,vg, ue,vn)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot
  real :: strcmp(9), xlat, xlong, ug, vg, ue, vn

  along = cspanf( xlong - strcmp(2), -180., 180.)
  if (xlat.gt.89.985) then
  !*  North polar meteorological orientation: "north" along prime meridian
    rot = - strcmp(1) * along + xlong - 180.
  elseif (xlat.lt.-89.985) then
  !*  South polar meteorological orientation: "north" along prime meridian
    rot = - strcmp(1) * along - xlong
  else
    rot = - strcmp(1) * along
  endif
  slong = sin( radpdg * rot )
  clong = cos( radpdg * rot )
  xpolg = slong * strcmp(5) + clong * strcmp(6)
  ypolg = clong * strcmp(5) - slong * strcmp(6)
  ue = ypolg * ug - xpolg * vg
  vn = ypolg * vg + xpolg * ug
  return
end subroutine cg2cll

subroutine cg2cxy (strcmp, x,y, ug,vg, ue,vn)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real :: strcmp(9) , x, y, ug, vg, ue, vn

  real :: clong, radial, rot, slong, xlat, xlong
  real(kind=dp) :: xpolg,ypolg,temp,xi0,eta0,xi,eta

  xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth
  eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth
  xi = xi0 * strcmp(5) - eta0 * strcmp(6)
  eta = eta0 * strcmp(5) + xi0 * strcmp(6)
  radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta)
  if (radial.gt.strcmp(8)) then
  !*  Case north of 89 degrees. Meteorological wind direction definition
  !*      changes.
    call cnxyll(strcmp, xi,eta, xlat,xlong)
  !*  North polar meteorological orientation: "north" along prime meridian
    rot = strcmp(1) * (xlong - strcmp(2)) - xlong - 180.
    slong = - sin( radpdg * rot )
    clong = cos( radpdg * rot )
    xpolg = slong * strcmp(5) + clong * strcmp(6)
    ypolg = clong * strcmp(5) - slong * strcmp(6)
  else if (radial.lt.strcmp(9)) then
  !*  Case south of -89 degrees. Meteorological wind direction definition
  !*      changes.
    call cnxyll(strcmp, xi,eta, xlat,xlong)
  !*  South polar meteorological orientation: "north" along prime meridian
    rot = strcmp(1) * (xlong - strcmp(2)) + xlong
    slong = - sin( radpdg * rot )
    clong = cos( radpdg * rot )
    xpolg = slong * strcmp(5) + clong * strcmp(6)
    ypolg = clong * strcmp(5) - slong * strcmp(6)
  else
  !* Normal case. Meteorological direction of wind related to true north.
    xpolg = strcmp(6) - strcmp(1) * xi0
    ypolg = strcmp(5) - strcmp(1) * eta0
    temp = sqrt ( xpolg ** 2 + ypolg ** 2 )
    xpolg = xpolg / temp
    ypolg = ypolg / temp
  end if
  ue = ( ypolg * ug - xpolg * vg )
  vn = ( ypolg * vg + xpolg * ug )
  return
end subroutine cg2cxy

real function cgszll (strcmp, xlat,xlong)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real :: strcmp(9), xlat, xlong

  real(kind=dp) :: slat,ymerc,efact

  if (xlat .gt. 89.985) then
  !* Close to north pole
    if (strcmp(1) .gt. 0.9999) then
  !* and to gamma == 1.
      cgszll = 2. * strcmp(7)
      return
    endif
    efact = cos(radpdg * xlat)
    if (efact .le. 0.) then
      cgszll = 0.
      return
    else
      ymerc = - log( efact /(1. + sin(radpdg * xlat)))
    endif
  else if (xlat .lt. -89.985) then
  !* Close to south pole
    if (strcmp(1) .lt. -0.9999) then
  !* and to gamma == -1.0
      cgszll = 2. * strcmp(7)
      return
    endif
    efact = cos(radpdg * xlat)
    if (efact .le. 0.) then
      cgszll = 0.
      return
    else
      ymerc = log( efact /(1. - sin(radpdg * xlat)))
    endif
  else
  slat = sin(radpdg * xlat)
  ymerc = log((1. + slat) / (1. - slat))/2.
  !efact = exp(ymerc)
  !cgszll = 2. * strcmp(7) * exp (strcmp(1) * ymerc)
  !c			 / (efact + 1./efact)
  endif
  cgszll = strcmp(7) * cos(radpdg * xlat) * exp(strcmp(1) *ymerc)
  return
end function cgszll

real function cgszxy (strcmp, x,y)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real :: strcmp(9) , x, y
  real(kind=dp) :: ymerc,efact, radial, temp
  real(kind=dp) :: xi0,eta0,xi,eta


  xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth
  eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth
  xi = xi0 * strcmp(5) - eta0 * strcmp(6)
  eta = eta0 * strcmp(5) + xi0 * strcmp(6)
  radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta)
  efact = strcmp(1) * radial
  if (efact .gt. almst1) then
    if (strcmp(1).gt.almst1) then
      cgszxy = 2. * strcmp(7)
    else
      cgszxy = 0.
    endif
    return
  endif
  if (abs(efact) .lt. 1.e-2) then
    temp = (efact / (2. - efact) )**2
    ymerc = radial / (2. - efact) * (1.    + temp * &
         (1./3. + temp * &
         (1./5. + temp * &
         (1./7. ))))
  else
    ymerc = - log( 1. - efact ) /2. /strcmp(1)
  endif
  if (ymerc .gt. 6.) then
    if (strcmp(1) .gt. almst1) then
      cgszxy = 2. * strcmp(7)
    else
      cgszxy = 0.
    endif
  else if (ymerc .lt. -6.) then
    if (strcmp(1) .lt. -almst1) then
      cgszxy = 2. * strcmp(7)
    else
      cgszxy = 0.
    endif
  else
    efact = exp(ymerc)
    cgszxy = 2. * strcmp(7) * exp (strcmp(1) * ymerc) &
         / (efact + 1./efact)
  endif
  return
end function cgszxy

subroutine cll2xy (strcmp, xlat,xlong, x,y)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  implicit none

  real :: strcmp(9) , xlat, xlong, x, y, xi, eta

  call cnllxy(strcmp, xlat,xlong, xi,eta)
  x = strcmp(3) + rearth/strcmp(7) * &
       (xi * strcmp(5) + eta * strcmp(6) )
  y = strcmp(4) + rearth/strcmp(7) * &
       (eta * strcmp(5) - xi * strcmp(6) )
  return
end subroutine cll2xy

subroutine cnllxy (strcmp, xlat,xlong, xi,eta)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL
  !  main transformation routine from latitude-longitude to
  !  canonical (equator-centered, radian unit) coordinates

  use par_mod, only: dp

  implicit none

  real :: strcmp(9), xlat, xlong, xi, eta, &
       gdlong, sndgam, csdgam, rhog1
  real(kind=dp) :: gamma
  real(kind=dp) :: dlong,dlat,slat,mercy,gmercy

  gamma = strcmp(1)
  dlat = xlat
  dlong = cspanf(xlong - strcmp(2), -180., 180.)
  dlong = dlong * radpdg
  gdlong = gamma * dlong
  if (abs(gdlong) .lt. .01) then
  !  Code for gamma small or zero.  This avoids round-off error or divide-
  !  by zero in the case of mercator or near-mercator projections.
    gdlong = gdlong * gdlong
    sndgam = dlong * (1. - 1./6. * gdlong * &
         (1. - 1./20. * gdlong * &
         (1. - 1./42. * gdlong )))
    csdgam = dlong * dlong * .5 * &
         (1. - 1./12. * gdlong * &
         (1. - 1./30. * gdlong * &
         (1. - 1./56. * gdlong )))
  else
  ! Code for moderate values of gamma
    sndgam = sin (gdlong) /gamma
    csdgam = (1. - cos(gdlong) )/gamma /gamma
  endif
  slat = sin(radpdg * dlat)
  if ((slat .ge. almst1) .or. (slat .le. -almst1)) then
    eta = 1./strcmp(1)
    xi = 0.
    return
  endif
  mercy = .5 * log( (1. + slat) / (1. - slat) )
  gmercy = gamma * mercy
  if (abs(gmercy) .lt. .001) then
  !  Code for gamma small or zero.  This avoids round-off error or divide-
  !  by zero in the case of mercator or near-mercator projections.
    rhog1 = mercy * (1. - .5 * gmercy * &
         (1. - 1./3. * gmercy * &
         (1. - 1./4. * gmercy ) ) )
  else
  ! Code for moderate values of gamma
    rhog1 = (1. - exp(-gmercy)) / gamma
  endif
  eta = rhog1 + (1. - gamma * rhog1) * gamma * csdgam
  xi = (1. - gamma * rhog1 ) * sndgam
end subroutine cnllxy

subroutine cnxyll (strcmp, xi,eta, xlat,xlong)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL
  !  main transformation routine from canonical (equator-centered,
  !  radian unit) coordinates

  use par_mod, only: dp

  implicit none

  real :: strcmp(9), xlat, xlong, odist
  real(kind=dp) :: gamma,temp,arg1,arg2,ymerc,along,gxi,cgeta
  real(kind=dp) :: xi,eta

  gamma = strcmp(1)
  !  Calculate equivalent mercator coordinate
  odist = xi*xi + eta*eta
  arg2 = 2. * eta - gamma * (xi*xi + eta*eta)
  arg1 = gamma * arg2
  ! Change by A. Stohl to avoid problems close to the poles
  ! if (arg1 .ge. almst1) then
  !  distance to north (or south) pole is zero (or imaginary ;) )
  ! xlat = sign(90.,strcmp(1))
  ! xlong = strcmp(2)
  ! return
  ! endif
  if (abs(arg1) .lt. .01) then
  !  Code for gamma small or zero.  This avoids round-off error or divide-
  !  by zero in the case of mercator or near-mercator projections.
    temp = (arg1 / (2. - arg1) )**2
    ymerc = arg2 / (2. - arg1) * (1.    + temp * &
         (1./3. + temp * &
         (1./5. + temp * &
         (1./7. ))))
  else
  ! Code for moderate values of gamma
    ymerc = - log ( 1. - arg1 ) /2. / gamma
  endif
  ! Convert ymerc to latitude
  temp = exp( - abs(ymerc) )
  xlat = sign(atan2((1. - temp) * (1. + temp), 2. * temp), ymerc)
  ! Find longitudes
  gxi = gamma*xi
  cgeta = 1. - gamma * eta
  if ( abs(gxi) .lt. .01*cgeta ) then
  !  Code for gamma small or zero.  This avoids round-off error or divide-
  !  by zero in the case of mercator or near-mercator projections.
    temp = ( gxi /cgeta )**2
    along = xi / cgeta * (1.    - temp * &
         (1./3. - temp * &
         (1./5. - temp * &
         (1./7.   ))))
  else
  ! Code for moderate values of gamma
    along = atan2( gxi, cgeta) / gamma
  endif
  xlong = sngl(strcmp(2) + dgprad * along)
  xlat = xlat * dgprad
  return
end subroutine cnxyll

subroutine cpolll (strcmp, xlat,xlong, enx,eny,enz)
  !*  Written on 11/23/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot
  real :: strcmp(9), xlat, xlong, enx, eny, enz, clat

  along = cspanf( xlong - strcmp(2), -180., 180.)
  rot = - strcmp(1) * along
  slong = sin( radpdg * rot )
  clong = cos( radpdg * rot )
  xpolg = slong * strcmp(5) + clong * strcmp(6)
  ypolg = clong * strcmp(5) - slong * strcmp(6)
  clat = cos(radpdg * xlat)
  enx = clat * xpolg
  eny = clat * ypolg
  enz = sin(radpdg * xlat)
  return
end subroutine cpolll

subroutine cpolxy (strcmp, x,y, enx,eny,enz)
  !*  Written on 11/26/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real :: strcmp(9) , x, y, enx, eny, enz
  real(kind=dp) :: xpol,ypol,temp,xi0,eta0,xi,eta,radial
  real(kind=dp) :: temp2,ymerc,arg,oarg,clat

  xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth
  eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth
  xi = xi0 * strcmp(5) - eta0 * strcmp(6)
  eta = eta0 * strcmp(5) + xi0 * strcmp(6)
  radial = 2. * eta -  strcmp(1) * (xi*xi + eta*eta)
  temp = strcmp(1) * radial
  if (temp .ge. 1.) then
    enx = 0.
    eny = 0.
    enz = sign(1.,strcmp(1))
    return
  endif
  if (abs(temp).lt.1.e-2) then
    temp2 = (temp / (2. - temp))**2
    ymerc = radial / (2. - temp) * (1. + temp2 * &
         (1./3. + temp2 * &
         (1./5. + temp2 * &
         (1./7.))))
  else
    ymerc = -.5 * log(1. - temp) / strcmp(1)
  endif
  arg = exp( ymerc )
  oarg = 1./arg
  clat = 2./(arg + oarg)
  enz = (arg - oarg) * clat /2.
  temp = clat / sqrt(1. - temp)
  xpol = - xi * strcmp(1) * temp
  ypol = (1. - eta * strcmp(1) ) * temp
  enx = xpol * strcmp(5) + ypol * strcmp(6)
  eny = ypol * strcmp(5) - xpol * strcmp(6)
  return
end subroutine cpolxy

real function cspanf (value, begin, end)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL
  !* real function cspanf returns a value in the interval (begin,end]
  !* which is equivalent to value, mod (end - begin).  It is used to
  !* reduce periodic variables to a standard range.  It adjusts for the
  !* behavior of the mod function which provides positive results for
  !* positive input, and negative results for negative input
  !* input:
  !*       value - real number to be reduced to the span
  !*       begin - first value of the span
  !*       end   - last value of the span
  !* returns:
  !*       the reduced value
  !* examples:
  !*      along = cspanf(xlong, -180., +180.)
  !*      dir  = cspanf(angle, 0., 360.)

  implicit none

  real :: first,last, value, begin, end, val

  first = min(begin,end)
  last = max(begin,end)
  val = mod( value - first , last - first)
  if ( val .le. 0.) then
    cspanf = val + last
  else
    cspanf = val + first
  endif
  return
end function cspanf

subroutine cxy2ll (strcmp, x,y, xlat,xlong)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  use par_mod, only: dp

  implicit none

  real(kind=dp) :: xi0,eta0,xi,eta
  real :: strcmp(9), x, y, xlat, xlong

  xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth
  eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth
  xi = xi0 * strcmp(5) - eta0 * strcmp(6)
  eta = eta0 * strcmp(5) + xi0 * strcmp(6)
  call cnxyll(strcmp, xi,eta, xlat,xlong)
  xlong = cspanf(xlong, -180., 180.)
  return
end subroutine cxy2ll

real function eqvlat (xlat1,xlat2)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  implicit none

  real :: xlat1, xlat2, x, ssind, sinl1, sinl2, al1, al2, tau

  ssind(x) = sin (radpdg*x)
  sinl1 = ssind (xlat1)
  sinl2 = ssind (xlat2)
  if (abs(sinl1 - sinl2) .gt. .001) then
    al1 = log((1. - sinl1)/(1. - sinl2))
    al2 = log((1. + sinl1)/(1. + sinl2))
  else
  !  Case lat1 near or equal to lat2
    tau = - (sinl1 - sinl2)/(2. - sinl1 - sinl2)
    tau = tau*tau
    al1  = 2. / (2. - sinl1 - sinl2) * (1.    + tau * &
         (1./3. + tau * &
         (1./5. + tau * &
         (1./7.))))
    tau =   (sinl1 - sinl2)/(2. + sinl1 + sinl2)
    tau = tau*tau
    al2  = -2. / (2. + sinl1 + sinl2) * (1.    + tau * &
         (1./3. + tau * &
         (1./5. + tau * &
         (1./7.))))
  endif
  eqvlat = asin((al1 + al2) / (al1 - al2))/radpdg
  return
end function eqvlat

subroutine stcm1p(strcmp, x1,y1, xlat1,xlong1, &
       xlatg,xlongg, gridsz, orient)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  implicit none

  integer :: k
  real :: strcmp(9), x1, y1, xlat1, xlong1, turn, orient, &
       xlatg, xlongg, gridsz, x1a, y1a

  do k=3,4
    strcmp (k) = 0.
  enddo
    turn = radpdg * (orient - strcmp(1) * &
         cspanf(xlongg - strcmp(2), -180., 180.) )
  strcmp (5) = cos (turn)
  strcmp (6) = - sin (turn)
  strcmp (7) = 1.
  strcmp (7) = gridsz * strcmp(7) &
       / cgszll(strcmp, xlatg, strcmp(2))
  call cll2xy (strcmp, xlat1,xlong1, x1a,y1a)
  strcmp(3) = strcmp(3) + x1 - x1a
  strcmp(4) = strcmp(4) + y1 - y1a
  return
end subroutine stcm1p

subroutine stcm2p(strcmp, x1,y1, xlat1,xlong1, &
       x2,y2, xlat2,xlong2)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  implicit none

  real :: strcmp(9), x1, y1, xlat1, xlong1, &
       x2, y2, xlat2, xlong2

  integer :: k
  real    :: x1a, y1a, x2a, y2a, den, dena

  do k=3,6
    strcmp (k) = 0.
  enddo
  strcmp (5) = 1.
  strcmp (7) = 1.
  call cll2xy (strcmp, xlat1,xlong1, x1a,y1a)
  call cll2xy (strcmp, xlat2,xlong2, x2a,y2a)
  den = sqrt( (x1 - x2)**2 + (y1 - y2)**2 )
  dena = sqrt( (x1a - x2a)**2 + (y1a - y2a)**2 )
  strcmp(5) = ((x1a - x2a)*(x1 - x2) + (y1a - y2a) * (y1 - y2)) &
       /den /dena
  strcmp(6) = ((y1a - y2a)*(x1 - x2) - (x1a - x2a) * (y1 - y2)) &
       /den /dena
  strcmp (7) = strcmp(7) * dena / den
  call cll2xy (strcmp, xlat1,xlong1, x1a,y1a)
  strcmp(3) = strcmp(3) + x1 - x1a
  strcmp(4) = strcmp(4) + y1 - y1a
  return
end subroutine stcm2p

!*  General conformal map routines for meteorological modelers
!*  written on 3/31/94 by

!* Dr. Albion Taylor
!* NOAA / OAR / ARL                  phone: (301) 713-0295 x 132
!* rm. 3151, 1315 east-west highway  fax:   (301) 713-0119
!* silver spring, md 20910           e-mail: adtaylor@arlrisc.ssmc.noaa.gov

!*  subroutine stlmbr (strcmp, tnglat, clong)
!*    This routine initializes the map structure array strcmp to
!*    the form of a specific map projection
!*  inputs:
!*    tnglat - the latitude at which the projection will be tangent
!*             to the earth.  +90. For north polar stereographic,
!*             -90. for south polar stereographic, 0. For mercator,
!*             and other values for lambert conformal.
!*             -90 <= tnglat <= 90.
!*    clong -  a longitude in the region under consideration.  Longitudes
!*             between clong-180. and clong+180.  Will be mapped in one
!*             connected region
!*  outputs:
!*    strcmp - a 9-value map structure array for use with subsequent
!*             calls to the coordinate transform routines.
!*
!*  real function eqvlat (xlat1,xlat2)
!*    This function is provided to assist in finding the tangent latitude
!*    equivalent to the 2-reference latitude specification in the legend
!*    of most lambert conformal maps.  If the map specifies "scale
!*    1:xxxxx true at 40n and 60n", then eqvlat(40.,60.) will return the
!*    equivalent tangent latitude.
!*  inputs:
!*    xlat1,xlat2:  the two latitudes specified in the map legend
!*  returns:
!*    the equivalent tangent latitude
!*  example:  call stlmbr(strcmp, eqvlat(40.,60.), 90.)

!*  subroutine stcm2p (strcmp, x1,y1, xlat1,xlong1,
!*          x2,y2, xlat2,xlong2)
!*  subroutine stcm1p (strcmp, x1,y1, xlat1,xlong1,
!*          xlatg,xlongg, gridsz, orient)
!*    These routines complete the specification of the map structure
!*    array by conforming the map coordinates to the specifications
!*    of a particular grid.  Either stcm1p or stcm2p must be called,
!*    but not both
!*  inputs:
!*    strcmp - a 9-value map structure array, set to a particular map
!*             form by a previous call to stlmbr
!*    for stcm2p:
!*      x1,y1, x2,y2 - the map coordinates of two points on the grid
!*      xlat1,xlong1, xlat2,xlong2 - the geographic coordinates of the
!*             same two points
!*    for stcm1p:
!*      x1,y1 - the map coordinates of one point on the grid
!*      xlat1,xlong1 - the geographic coordinates of the same point
!*      xlatg,xlongg - latitude and longitude of reference point for
!*             gridsz and orientation specification.
!*      gridsz - the desired grid size in kilometers, at xlatg,xlongg
!*      orient - the angle, with respect to north, of a y-grid line, at
!*             the point xlatg,xlongg
!*  outputs:
!*    strcmp - a 9-value map structure array, fully set for use by
!*             other subroutines in this system

!*  subroutine cll2xy (strcmp, xlat,xlong, x,y)
!*  subroutine cxy2ll (strcmp, x,y, xlat,xlong)
!*     these routines convert between map coordinates x,y
!*     and geographic coordinates xlat,xlong
!*  inputs:
!*     strcmp(9) - 9-value map structure array
!*     for cll2xy:  xlat,xlong - geographic coordinates
!*     for cxy2ll:  x,y - map coordinates
!*  outputs:
!*     for cll2xy:  x,y - map coordinates
!*     for cxy2ll:  xlat,xlong - geographic coordinates

!*  subroutine cc2gxy (strcmp, x,y, ue,vn, ug,vg)
!*  subroutine cg2cxy (strcmp, x,y, ug,vg, ue,vn)
!*  subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg)
!*  subroutine cg2cll (strcmp, xlat,xlong, ug,vg, ue,vn)
!*     These subroutines convert vector wind components from
!*     geographic, or compass, coordinates, to map or
!*     grid coordinates.  The site of the wind to be
!*     converted may be given either in geographic or
!*     map coordinates.  Wind components are all in kilometers
!*     per hour, whether geographic or map coordinates.
!*  inputs:
!*    strcmp(9) - 9-value map structure array
!*    for cc2gxy and cg2cxy:  x,y        -  map coordinates of site
!*    for cc2gll and cg2cll:  xlat,xlong -  geographic coordinates of site
!*    for cc2gxy and cc2gll:  ue,vn - east and north wind components
!*    for cg2cxy and cg2cll:  ug,vg - x- and y- direction wind components
!*  outputs:
!*    for cc2gxy and cc2gll:  ug,vg - x- and y- direction wind components
!*    for cg2cxy and cg2cll:  ue,vn - east and north wind components

!*  subroutine ccrvxy (strcmp, x, y,       gx,gy)
!*  subroutine ccrvll (strcmp, xlat,xlong, gx,gy)
!*    These subroutines return the curvature vector (gx,gy), as referenced
!*    to map coordinates, induced by the map transformation.  When
!*    non-linear terms in wind speed are important, a "geodesic" force
!*    should be included in the vector form [ (u,u) g - (u,g) u ] where the
!*    inner product (u,g) is defined as ux*gx + uy*gy.
!*  inputs:
!*    strcmp(9) - 9-value map structure array
!*    for ccrvxy:  x,y        -  map coordinates of site
!*    for ccrvll:  xlat,xlong -  geographic coordinates of site
!*  outputs:
!*    gx,gy       - vector coefficients of curvature, in units radians
!*                  per kilometer

!*  real function cgszll (strcmp, xlat,xlong)
!*  real function cgszxy (strcmp, x,y)
!*    These functions return the size, in kilometers, of each unit of
!*    motion in map coordinates (grid size).  The grid size at any
!*    location depends on that location; the position may be given in
!*    either map or geographic coordinates.
!*  inputs:
!*    strcmp(9) - 9-value map structure array
!*    for cgszxy:  x,y        -  map coordinates of site
!*    for cgszll:  xlat,xlong -  geographic coordinates of site
!*  returns:
!*    gridsize in kilometers at given site.

!*  subroutine cpolxy (strcmp, x,y, enx,eny,enz)
!*  subroutine cpolll (strcmp, xlat,xlong, enx,eny,enz)
!*    These subroutines provide 3-d vector components of a unit vector
!*    in the direction of the north polar axis.  When multiplied
!*    by twice the rotation rate of the earth (2 * pi/24 hr), the
!*    vertical component yields the coriolis factor.
!*  inputs:
!*    strcmp(9) - 9-value map structure array
!*    for cpolxy:  x,y        -  map coordinates of site
!*    for cpolll:  xlat,xlong -  geographic coordinates of site
!*  returns:
!*    enx,eny,enz the direction cosines of a unit vector in the
!*    direction of the rotation axis of the earth

!*  subroutine cnllxy (strcmp, xlat,xlong, xi,eta)
!*  subroutine cnxyll (strcmp, xi,eta, xlat,xlong)
!*    These subroutines perform the underlying transformations from
!*    geographic coordinates to and from canonical (equator centered)
!*    coordinates.  They are called by cxy2ll and cll2xy, but are not
!*    intended to be called directly

!*  real function cspanf (value, begin, end)
!*    This function assists other routines in providing a longitude in
!*    the proper range.  It adds to value whatever multiple of
!*    (end - begin) is needed to return a number begin < cspanf <= end

subroutine stlmbr(strcmp, tnglat, xlong)
  !*  Written on 3/31/94 by Dr. Albion Taylor  NOAA / OAR / ARL

  implicit none

  real :: strcmp(9), tnglat, xlong

  real :: eta, xi

  strcmp(1) = sin(radpdg * tnglat)
  !*  gamma = sine of the tangent latitude
  strcmp(2) = cspanf( xlong, -180., +180.)
  !* lambda_0 = reference longitude
  strcmp(3) = 0.
  !* x_0 = x- grid coordinate of origin (xi,eta) = (0.,0.)
  strcmp(4) = 0.
  !* y_0 = y-grid coordinate of origin (xi,eta) = (0.,0.)
  strcmp(5) = 1.
  !* Cosine of rotation angle from xi,eta to x,y
  strcmp(6) = 0.
  !* Sine of rotation angle from xi,eta to x,y
  strcmp(7) = rearth
  !* Gridsize in kilometers at the equator
  call cnllxy(strcmp, 89.,xlong, xi,eta)
  strcmp(8) = 2. * eta - strcmp(1) * eta * eta
  !* Radial coordinate for 1 degree from north pole
  call cnllxy(strcmp, -89.,xlong, xi,eta)
  strcmp(9) = 2. * eta - strcmp(1) * eta * eta
  !* Radial coordinate for 1 degree from south pole
  return
end subroutine stlmbr

end module cmapf_mod