pro eqtoll, phi_eq,lambda_eq,phi,lambda,phi_pole,lambda_pole         
;+  
;  Procedure eqtoll ------------------------------------------------- 
;                                                                     
;  Purpose:  Calculates latitude and longitude on standard grid       
;            from input arrays of latitude and longitude on           
;            equatorial latitude-longitude (eq) grid used             
;            in regional models. Both input and output latitudes      
;            and longitudes are in degrees.                           
;                                                                     
;  Written by Pete Clark, translating A. Dickinson FORTRAN to WAVE    
;  24/09/2002 - Name of routine sgn changed to apply_sgn                                       
;                                                                     
;  Documentation: The transformation formulae are described in        
;                 unified model on-line documentation paper S1.       
;                                                                     
; Logical components covered : S131                                   
;                                                                     
; Project task :                                                      
;                                                                     
; External documentation:                                             
;                                                                     
;  Arguments:-------------------------------------------------------- 
;
; phi_eq            !in  latitude in equatorial lat-lon coords  model    
; lambda_eq         !in  longitude in equatorial lat-lon coords    
; phi               !out latitude                               normal
; lambda            !out longitude (0 =< lon < 360)                
; phi_pole          !in  latitude of equatorial lat-lon pole        
; lambda_pole       !in  longitude of equatorial lat-lon pole    
;-  

 
  small=1.0e-7                                           
                                                                        
  pi=!pi     ; pi                                    
  pi_over_180 =!dtor       ; conversion factor degrees to radians  
  recip_pi_over_180 = !radeg; conversion factor radians to degrees       
                                                                        
; 1. initialise local constants                                        
; latitude of zeroth meridian                                           
  if lambda_pole gt 90 or lambda_pole lt -90 then begin
    lambda_zero=lambda_pole+180. 
  endif else begin                                  
    lambda_zero=lambda_pole
  endelse                                   
; 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) 
  
; print, 'spole', sin_phi_pole, cos_phi_pole                         
                                                                        
; 2. transform from equatorial to standard latitude-longitude          
                                                                        
                                                                        
; scale eq longitude to range -180 to +180 degs                         
                                                                        
  e_lambda=lambda_eq
                                             
  i=where(e_lambda gt 180.0,count) 
  if count gt 0 then e_lambda(i)=e_lambda(i)-360.0                    
  i=where(e_lambda lt -180.0,count) 
  if count gt 0 then e_lambda(i)=e_lambda(i)+360.0    
;  info,e_lambda                
                                                                        
; convert eq latitude & longitude to radians                            
                                                                        
  e_lambda=pi_over_180*e_lambda                                     
;  info,e_lambda                
  e_phi=pi_over_180*phi_eq                                       
;  info,e_phi              
                                                                        
; compute latitude using equation (4.7)                                 
                                                                        
  arg=cos_phi_pole*cos(e_lambda)*cos(e_phi) + sin(e_phi)*sin_phi_pole  
;  info,arg 
  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                                                 
;  info,arg                     
  a_phi=asin(arg)                                                   
;  info,a_phi                                 
  phi=recip_pi_over_180*a_phi   
;  info,phi                                 
                                                                        
; compute longitude using equation (4.8)                                
                                                                        
  term1 =(cos(e_phi)*cos(e_lambda)*sin_phi_pole - sin(e_phi)*cos_phi_pole)                                  
  term2=cos(a_phi) 
;  info,term1,term2                                                 
  i1=where(term2 lt small,count1)                                         
  i2=where(term2 ge small,count2)  
;  info,count1,count2
  a_lambda=term1 ;  dummy to set up space
  
  if count1 gt 0 then begin                                     
    a_lambda(i1)=0.0  
    info,a_lambda
  endif
                                                
  if count2 gt 0 then begin  
;    print,'Non zero'                                   
    arg=term1/term2 
;    info,arg
;    print,term1-term2
    i= where(arg ge 1.0-small,count)  
                                                 
    if count gt 0 then arg(i) = 1.0 -small  
                                              
    a_lambda(i2)=recip_pi_over_180*acos(arg(i2)) 
    if count gt 0 then a_lambda(i)=0
;    info,a_lambda
    a_lambda(i2) = apply_sgn(a_lambda(i2), e_lambda(i2))  
                                  
    a_lambda(i2) = a_lambda(i2)+lambda_zero(i2)                                   
;    info,a_lambda
  endif
                                                                        
; scale longitude to range 0 to 360 degs                                
                                                                        
  i1=where(a_lambda ge 360.0,count) 
  if count gt 0 then a_lambda(i1)=a_lambda(i1)-360.0                     
  i1=where(a_lambda lt 0.0,count) 
  if count gt 0 then a_lambda(i1)=a_lambda(i1)+360.0                     
  lambda=a_lambda                                                
                                                                        
  return                                                            
end