m_oldtonew.f90 1.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102
  1. # 0 "<stdin>"
  2. # 0 "<built-in>"
  3. # 0 "<command-line>"
  4. # 1 "/usr/include/stdc-predef.h" 1 3 4
  5. # 17 "/usr/include/stdc-predef.h" 3 4
  6. # 2 "<command-line>" 2
  7. # 1 "<stdin>"
  8. # 10 "<stdin>"
  9. module m_oldtonew
  10. use m_confmap
  11. implicit none
  12. contains
  13. ! this routine performes a conformal mapping of the old to the new
  14. ! coordinate system
  15. !
  16. subroutine oldtonew(lat_o, lon_o, lat_n, lon_n)
  17. real(8), intent(in) :: lat_o, lon_o
  18. real(8), intent(out) :: lat_n, lon_n
  19. real :: theta, phi, psi, mu
  20. complex :: z, w
  21. if (.not. confmap_initialised) then
  22. print *, 'ERROR: oldtonew(): confmap not initialised'
  23. stop
  24. end if
  25. ! transform to spherical coordinates
  26. !
  27. theta = mod(lon_o * rad + 3.0 * pi_1, 2.0 * pi_1) - pi_1
  28. phi = pi_2 - lat_o * rad
  29. ! transform to the new coordinate system
  30. !
  31. if (abs(phi - pi_1) < epsil) then
  32. mu = mu_s
  33. psi = psi_s
  34. elseif (abs(phi - phi_b) < epsil .and. abs(theta - theta_b) < epsil) then
  35. mu = 0.0
  36. psi = pi_1
  37. else
  38. z = tan(0.5 * phi) * exp(imagone * theta)
  39. w = (z - ac) * cmnb / ((z - bc) * cmna)
  40. mu = atan2(aimag(w), real(w))
  41. psi = 2.0 * atan(abs(w))
  42. endif
  43. ! transform to lat/lon coordinates
  44. !
  45. lat_n = (pi_2 - psi) * deg
  46. lon_n = mu * deg
  47. end subroutine oldtonew
  48. end module m_oldtonew