;+  pro lltoeq-------------------------------------------------
;
;  written by a. dickinson 11/8/89
;
;  reviewed by r. swinbank 19/01/90
;
;  purpose:
;            calculates latitude and longitude on equatorial
;            latitude-longitude (eq) grid used in regional
;            models from input arrays of latitude and
;            longitude on standard grid. both input and output
;            latitudes and longitudes are in degrees.
;
;  documentation:
;            the transformation formulae are described in
;            unified model on-line documentation paper s1.
;  24/09/2002 - Name of routine sgn changed to apply_sgn                                       
;
;  ------------------------------------------------------------------
;  Arguments:-------------------------------------------------------- 
;
; phi               !in latitude                                  
; lambda            !in longitude (0 =< lon < 360)                
; lambda_eq         !out longitude in equatorial lat-lon coords    
; phi_eq            !out latitude in equatorial lat-lon coords     
; phi_pole          !in  latitude of equatorial lat-lon pole       
; lambda_pole       !in  longitude of equatorial lat-lon pole    
;-  
pro lltoeq, phi,lambda,phi_eq,lambda_eq,phi_pole,lambda_pole

  small=1.0e-6
  pi=3.14159265358979323846     ; pi                                    
  pi_over_180 =pi/180.0        ; conversion factor degrees to radians  
  recip_pi_over_180 = 180.0/pi ; conversion factor radians to degrees       
; latitude of zeroth meridian
  lambda_zero=lambda_pole+180.
; sine and cosine of latitude of eq pole
  sin_phi_pole=sin(pi_over_180*phi_pole)
  cos_phi_pole=cos(pi_over_180*phi_pole)

; transform from standard to equatorial latitude-longitude

; scale longitude to range -180 to +180 degs

  a_lambda=lambda-lambda_zero
  i=where(a_lambda gt 180.0,count) 
  if count gt 0 then a_lambda(i)=a_lambda(i)-360.0                    
  i=where(a_lambda le -180.0,count) 
  if count gt 0 then a_lambda(i)=a_lambda(i)+360.0                    

; convert latitude & longitude to radians

  a_lambda=pi_over_180*a_lambda
  a_phi=pi_over_180*phi

; compute eq latitude using equation (4.4)

  arg=-cos_phi_pole*cos(a_lambda)*cos(a_phi) + sin(a_phi)*sin_phi_pole
     
  i=where(arg gt 1.0, count)
  if count gt 0 then arg(i)=1.0
  i=where(arg lt -1.0, count)
  if count gt 0 then arg(i)=-1.0
  e_phi=asin(arg)
  phi_eq=recip_pi_over_180*e_phi

; compute eq longitude using equation (4.6)

  term1 =cos(a_phi)*cos(a_lambda)*sin_phi_pole + sin(a_phi)*cos_phi_pole
  term2 =cos(e_phi)
  i1=where(term2 lt small,count1)
  i2=where(term2 ge small,count2)
  e_lambda=lambda
  if count1 gt 0 then begin
    e_lambda(i1)=0.0
  endif
  if count2 gt 0 then begin
    arg=term1(i2)/term2(i2)
    i=where(arg gt 1.0, count)
    if count gt 0 then arg(i)=1.0
    i=where(arg lt -1.0, count)
    if count gt 0 then arg(i)=-1.0
    e_lambda(i2)=recip_pi_over_180*acos(arg)
    e_lambda(i2)=apply_sgn(e_lambda(i2),a_lambda(i2))
  endif

; scale longitude to range 0 to 360 degs

  i=where(e_lambda ge 360.0,count) 
  if count gt 0 then e_lambda(i)=e_lambda(i)-360.0                    
  i=where(e_lambda lt 0.0,count) 
  if count gt 0 then e_lambda(i)=e_lambda(i)+360.0                    
  lambda_eq=e_lambda
  return
end