Index: LMDZ5/branches/testing/libf/dynphy_lonlat/grid_atob_m.f90
===================================================================
--- LMDZ5/branches/testing/libf/dynphy_lonlat/grid_atob_m.f90	(revision 2532)
+++ LMDZ5/branches/testing/libf/dynphy_lonlat/grid_atob_m.f90	(revision 2533)
@@ -29,16 +29,14 @@
 ! Local variables:
   CHARACTER(LEN=256) :: modname="fine2coarse"
-  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1,m2, n1,n2, mx,my, nn, i,j
-  LOGICAL :: li, lo, first=.TRUE.
-  REAL    :: inc, cpa, spa, crlo(SIZE(x_i))
-  REAL, SAVE :: pi, hpi
+  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn
+  LOGICAL :: li, lo
+  REAL    :: inc
   INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
   LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
-  REAL,    DIMENSION(SIZE(x_i),SIZE(y_i)) :: dist
+  REAL,    DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist
   REAL,    DIMENSION(SIZE(x_o))           :: a, b
   REAL,    DIMENSION(SIZE(y_o))           :: c, d
   REAL,    PARAMETER :: thresh=1.E-5
 !-------------------------------------------------------------------------------
-  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
   mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
   m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
@@ -93,7 +91,5 @@
       IF(found(io,jo)) CYCLE
 !      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
-      crlo(:)=COS(x_o(io)-x_i(:))     !--- COS of points 1 and 2 angle
-      cpa=COS(y_o(jo)); spa=SIN(y_o(jo))
-      DO j=1,ni; dist(:,j)=ACOS(spa*SIN(y_i(j))+cpa*COS(y_i(j))*crlo(:)); END DO
+      CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:))
       nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
       inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
@@ -256,4 +252,39 @@
 !-------------------------------------------------------------------------------
 
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance)
+!
+!-------------------------------------------------------------------------------
+! Author:  Laurent Li (december 30th 1996).
+! Purpose: Compute min. distance (along big circle) between 2 points in radians.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (radians)
+  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians)
+  REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (radians)
+!-------------------------------------------------------------------------------
+! Local variables:
+  LOGICAL, SAVE :: first=.TRUE.
+  REAL,    SAVE :: pi, hpi
+  REAL    :: pa, pb, cpa, spa, crlo(SIZE(rlon))
+  INTEGER :: i, j
+!-------------------------------------------------------------------------------
+  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
+  crlo(:)=COS(rf_lon-rlon(:))         !--- COS of points 1 and 2 angle
+  pa=hpi-rf_lat                       !--- North Pole - Point 1 distance
+  cpa=COS(pa); spa=SIN(pa)
+  DO j=1,SIZE(rlat)
+    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
+    distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:))
+  END DO
+
+END SUBROUTINE dist_sphe
+!
+!-------------------------------------------------------------------------------
+
 END MODULE grid_atob_m
 !