| 1 | ! $Id: occ_point.f90 2019 2009-01-14 10:20:26Z frhl $
|
|---|
| 2 |
|
|---|
| 3 | !****f* Coordinates/occ_point
|
|---|
| 4 | !
|
|---|
| 5 | ! NAME
|
|---|
| 6 | ! occ_point - Determine the occultation point
|
|---|
| 7 | !
|
|---|
| 8 | ! SYNOPSIS
|
|---|
| 9 | ! call occ_point(r_leo, r_gns, lat, lon, r_coc, roc, azimuth, undulation,
|
|---|
| 10 | ! cfile, efile)
|
|---|
| 11 | !
|
|---|
| 12 | ! DESCRIPTION
|
|---|
| 13 | ! This subroutine calculates the lowest occultation perigree point projected
|
|---|
| 14 | ! to the Earth's surface
|
|---|
| 15 | !
|
|---|
| 16 | ! INPUTS
|
|---|
| 17 | ! r_leo cartesian LEO position vector (relative to ECF frame)
|
|---|
| 18 | ! r_gns cartesian GPS position vector (relative to ECF frame)
|
|---|
| 19 | ! cfile path to geoid potential coefficients file (optional)
|
|---|
| 20 | ! efile path to geoid potential corrections file (optional)
|
|---|
| 21 | !
|
|---|
| 22 | ! OUTPUT
|
|---|
| 23 | ! lat Occultation point latitude
|
|---|
| 24 | ! lon Occultation point longitude
|
|---|
| 25 | ! r_coc Cartesian centre of curvature vector for occ point (ECF)
|
|---|
| 26 | ! roc Radius of curvature value for occultation point
|
|---|
| 27 | ! azimuth GPS to LEO azimuth direction wrt true North (deg)
|
|---|
| 28 | ! undulation Difference between ellipsoid (WGS-84) and EGM-96 geoid (m)
|
|---|
| 29 | !
|
|---|
| 30 | ! SEE ALSO
|
|---|
| 31 | !
|
|---|
| 32 | ! REFERENCES
|
|---|
| 33 | !
|
|---|
| 34 | ! AUTHOR
|
|---|
| 35 | ! Met Office, Exeter, UK.
|
|---|
| 36 | ! Any comments on this software should be given via the ROM SAF
|
|---|
| 37 | ! Helpdesk at http://www.romsaf.org
|
|---|
| 38 | !
|
|---|
| 39 | ! COPYRIGHT
|
|---|
| 40 | ! (c) EUMETSAT. All rights reserved.
|
|---|
| 41 | ! For further details please refer to the file COPYRIGHT
|
|---|
| 42 | ! which you should have received as part of this distribution.
|
|---|
| 43 | !
|
|---|
| 44 | !****
|
|---|
| 45 |
|
|---|
| 46 | subroutine occ_point(r_leo, r_gns, lat, lon, r_coc, roc, azimuth, &
|
|---|
| 47 | undulation, cfile,efile)
|
|---|
| 48 |
|
|---|
| 49 | ! 1.1 Declarations
|
|---|
| 50 | ! ----------------
|
|---|
| 51 |
|
|---|
| 52 | use typesizes, only: wp => EightByteReal
|
|---|
| 53 | use coordinates, not_this => occ_point
|
|---|
| 54 | use EarthMod, only: datum_hmsl
|
|---|
| 55 |
|
|---|
| 56 | implicit none
|
|---|
| 57 |
|
|---|
| 58 | real(wp), dimension(:,:), intent(in) :: r_leo ! LEO position vector (ECF)
|
|---|
| 59 | real(wp), dimension(:,:), intent(in) :: r_gns ! GPS position vector (ECF)
|
|---|
| 60 | real(wp), intent(out) :: lat ! Occultation point latitude
|
|---|
| 61 | real(wp), intent(out) :: lon ! Occultation point longitude
|
|---|
| 62 | real(wp), dimension(size(r_leo,2)), intent(out) :: r_coc ! Centre curvature
|
|---|
| 63 | real(wp), intent(out) :: roc ! Radius curvature
|
|---|
| 64 | real(wp), intent(out) :: azimuth ! Azimuth (deg)
|
|---|
| 65 | real(wp), intent(out) :: undulation ! Undulation
|
|---|
| 66 | character(len=*), optional, intent(in) :: cfile ! Coefficient file path
|
|---|
| 67 | character(len=*), optional, intent(in) :: efile ! Corrections file path
|
|---|
| 68 |
|
|---|
| 69 | real(wp), dimension(size(r_leo,1),size(r_leo,2)) :: perigee
|
|---|
| 70 | real(wp) :: slta ! Straight line tangent ht
|
|---|
| 71 | real(wp) :: ro ! Length of r_leo
|
|---|
| 72 | real(wp) :: alpha ! Angle r_leo and perigee
|
|---|
| 73 | real(wp) :: theta ! Cross section azimuth
|
|---|
| 74 | real(wp), dimension(size(r_leo,1)) :: height_per
|
|---|
| 75 | real(wp), dimension(size(r_leo,1)) :: lat_per
|
|---|
| 76 | real(wp), dimension(size(r_leo,1)) :: lon_per
|
|---|
| 77 | real(wp), dimension(3), parameter :: pa = (/0,0,1/) ! Polar axis
|
|---|
| 78 |
|
|---|
| 79 | integer:: i, iocc
|
|---|
| 80 |
|
|---|
| 81 | do i=1,size(r_leo,1)
|
|---|
| 82 |
|
|---|
| 83 | ! 1.2 Determine ray tangent points
|
|---|
| 84 | ! --------------------------------
|
|---|
| 85 |
|
|---|
| 86 | slta = impact_parameter(r_leo(i,:), r_gns(i,:))
|
|---|
| 87 | ro = Sqrt(Dot_Product(r_leo(i,:), r_leo(i,:)))
|
|---|
| 88 | alpha = acos(slta/ro)
|
|---|
| 89 |
|
|---|
| 90 | perigee(i,:) = rotate(r_leo(i,:), vector_product(r_leo(i,:), r_gns(i,:)), alpha) * (slta/ro)
|
|---|
| 91 |
|
|---|
| 92 | enddo
|
|---|
| 93 |
|
|---|
| 94 | ! 1.3 Convert cartesian to geodetic points
|
|---|
| 95 | ! ----------------------------------------
|
|---|
| 96 |
|
|---|
| 97 | call cart2geod(perigee, lat_per, lon_per, height_per)
|
|---|
| 98 |
|
|---|
| 99 | ! 1.4 Find the lowest perigee
|
|---|
| 100 | ! ---------------------------
|
|---|
| 101 |
|
|---|
| 102 | iocc = Sum(MinLoc(Abs(height_per)))
|
|---|
| 103 |
|
|---|
| 104 | ! 1.5 Define occultation point latitude and longitude
|
|---|
| 105 | ! ---------------------------------------------------
|
|---|
| 106 |
|
|---|
| 107 | lat = lat_per(iocc)
|
|---|
| 108 | lon = lon_per(iocc)
|
|---|
| 109 |
|
|---|
| 110 | ! 1.6 Determine occultation point centre of curvature and radius
|
|---|
| 111 | ! --------------------------------------------------------------
|
|---|
| 112 |
|
|---|
| 113 | ! 1.6.1 Cross-section azimuth
|
|---|
| 114 |
|
|---|
| 115 | theta = vector_angle( vector_product(perigee(iocc,:), pa), &
|
|---|
| 116 | vector_product(r_gns(iocc,:), r_leo(iocc,:)), &
|
|---|
| 117 | -perigee(iocc,:) )
|
|---|
| 118 |
|
|---|
| 119 | azimuth = theta * rad2deg
|
|---|
| 120 |
|
|---|
| 121 | if ( azimuth < 0.0_wp ) azimuth = azimuth + 360.0_wp
|
|---|
| 122 |
|
|---|
| 123 | ! 1.6.2 Compute curvature
|
|---|
| 124 |
|
|---|
| 125 | call curvature(lat, lon, theta, r_coc, roc)
|
|---|
| 126 |
|
|---|
| 127 | ! 1.7 Find undulation - height difference between local ellipsoid and geoid
|
|---|
| 128 | ! -------------------------------------------------------------------------
|
|---|
| 129 |
|
|---|
| 130 | if (present(cfile) .and. present(efile)) then
|
|---|
| 131 | call datum_hmsl("WGS84", (/lat, lon, 0.0_wp/), undulation, cfile, efile)
|
|---|
| 132 | else
|
|---|
| 133 | call datum_hmsl("WGS84", (/lat, lon, 0.0_wp/), undulation)
|
|---|
| 134 | endif
|
|---|
| 135 | if (undulation > -999999.000) undulation = -1.0_wp * undulation
|
|---|
| 136 |
|
|---|
| 137 | end subroutine occ_point
|
|---|
| 138 |
|
|---|
| 139 |
|
|---|
| 140 |
|
|---|
