bipolar_corner_latlon Subroutine

   pure subroutine bipolar_corner_latlon(lam_deg, s, phi_join, lon_pole, lat, lon)
      !! Map a logical cap location to geographic (lat, lon) in degrees.
      !!   `lam_deg` : pseudo-longitude (geographic lon at the join ring),
      !!               any real (wrapped internally); the i-direction.
      !!   `s`       : cap-row fraction, 0 = join ring (`lat = phi_join`),
      !!               1 = fold line; the j-direction.
      !!   `phi_join`: join latitude (deg); cap covers lat > phi_join.
      !!   `lon_pole`: longitude (deg) of the first cap pole; partner at
      !!               `lon_pole + 180`.
      real(wp), intent(in) :: lam_deg, s, phi_join, lon_pole
      real(wp), intent(out) :: lat, lon

      real(wp) :: a_focus, half, t, mu, xi, sgn
      real(wp) :: wr, wi, denr, deni, dnorm, zpr, zpi
      real(wp) :: lpr, zr, zi, rmag

      a_focus = tan(DEG2RAD*(90.0_wp - phi_join)*0.5_wp)

      ! mu from the join-continuity relation: lon_join = lon_pole + 2 atan(exp(mu)).
      half = DEG2RAD*(lam_deg - lon_pole)*0.5_wp
      t = tan(half)
      if (t > 0.0_wp) then
         sgn = 1.0_wp
         mu = log(t)
      else if (t < 0.0_wp) then
         sgn = -1.0_wp
         mu = log(-t)
      else
         ! lam == lon_pole (mod 360): a cap pole.  mu -> -inf; clamp finite
         ! so the corner sits arbitrarily close to the pole without an Inf.
         sgn = 1.0_wp
         mu = -50.0_wp
      end if

      ! xi rotates from the join (sgn*pi/2) to the fold (sgn*pi) as s: 0 -> 1.
      xi = sgn*(PI*0.5_wp + s*PI*0.5_wp)

      ! w = exp(mu) * exp(i*xi)
      wr = exp(mu)*cos(xi)
      wi = exp(mu)*sin(xi)

      ! Zp = A * (1 + w)/(1 - w)
      denr = 1.0_wp - wr
      deni = -wi
      dnorm = denr*denr + deni*deni
      if (dnorm <= 0.0_wp) then
         ! w == 1 (mu -> +inf): the other pole.  Place at phi_join on the
         ! partner meridian.
         lat = phi_join
         lon = lon_pole + 180.0_wp
         if (lon >= 180.0_wp) lon = lon - 360.0_wp
         return
      end if
      ! (1+w)/(1-w)
      zpr = ((1.0_wp + wr)*denr + wi*deni)/dnorm
      zpi = (wi*denr - (1.0_wp + wr)*deni)/dnorm
      zpr = a_focus*zpr
      zpi = a_focus*zpi

      ! Z = Zp * exp(i*lon_pole)
      lpr = DEG2RAD*lon_pole
      zr = zpr*cos(lpr) - zpi*sin(lpr)
      zi = zpr*sin(lpr) + zpi*cos(lpr)

      ! stereographic inverse
      rmag = sqrt(zr*zr + zi*zi)
      lat = 90.0_wp - RAD2DEG*2.0_wp*atan(rmag)
      if (zr == 0.0_wp .and. zi == 0.0_wp) then
         lon = lon_pole + 180.0_wp
      else
         lon = RAD2DEG*atan2(zi, zr)
      end if
      ! wrap lon to [-180, 180)
      lon = modulo(lon + 180.0_wp, 360.0_wp) - 180.0_wp
   end subroutine bipolar_corner_latlon