Index: /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90
===================================================================
--- /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90	(revision 2857)
+++ /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90	(revision 2858)
@@ -58,5 +58,4 @@
         zdistante  , zdistante2 , zdistnew  , zdistnew2  ,  &!         "
         zdeltai    , zdeltaj                                !         "
-!     REAL(wp),DIMENSION(:,:),ALLOCATABLE,SAVE :: zmask
      LOGICAL :: & 
         ll_overlap_sec_left = .FALSE. , ll_overlap_sec_right = .FALSE. ,&! temporary logical
@@ -75,5 +74,5 @@
         SouthPoint , NorthPoint, EstPoint , WestPoint  !        "
      !!---------------------------------------------------------------------
-     IF ( jsec==1 )THEN
+     IF( jsec==1 )THEN
           PRINT*,'                '
           PRINT*,'COMPUTE SECTIONS'
@@ -104,38 +103,18 @@
      SouthPoint  = POINT_SECTION( -1, -1 ) ; NorthPoint  = POINT_SECTION( -1, -1 )
      EstPoint    = POINT_SECTION( -1, -1 ) ; WestPoint   = POINT_SECTION( -1, -1 )
-   
-     !PRINT*,"min max tmask ",MINVAL(tmask),MAXVAL(tmask)
-
-     IF (jsec == 1)THEN
-         !Found the taller cell of the mesh in ocean (used in compsec)
-         zdistmesh=0.
-         DO jj=2,nlej
-            DO ji=2,nlei
-               IF( zdistmesh .LT. &
-                 !(e1t(ji,jj)/(cos(gphit(ji,jj))*110000)) *tmask(ji,jj,1 ) ) THEN
-                 !zdistmesh = e1t(ji,jj)/(cos(gphit(ji,jj))*110000)*tmask(ji,jj,1 )
+
+     IF( jsec == 1 )THEN
+        !Found the taller cell of the mesh in ocean (used in compsec)
+        zdistmesh=0.
+        DO jj=2,jpj
+           DO ji=2,jpi
+              IF( zdistmesh .LT. &
                  (e1t(ji,jj)/(cos(gphit(ji,jj))*110000))  ) THEN
                  zdistmesh = e1t(ji,jj)/(cos(gphit(ji,jj))*110000)
-               ENDIF
-            ENDDO
-         ENDDO
+              ENDIF
+           ENDDO
+        ENDDO
     
-         !Compute array zmask used later to avoid array overflow  
-         !ALLOCATE(zmask(jpi,jpj))
-         !zmask(:,:)=tmask(:,:,1)
-         !DO jj=3,nlcj-2
-         !   DO ji=3,nlci-2
-         !      zmask(ji,jj) =  MIN( 1.,     tmask(ji  ,jj  ,1)  &
-         !           + tmask(ji-1,jj  ,1) + tmask(ji+1,jj  ,1)  &
-         !           + tmask(ji-2,jj  ,1) + tmask(ji+2,jj  ,1)  &
-         !           + tmask(ji  ,jj-1,1) + tmask(ji  ,jj+1,1)  & 
-         !           + tmask(ji  ,jj-2,1) + tmask(ji  ,jj+2,1)  &
-         !           + tmask(ji-1,jj-1,1) + tmask(ji+1,jj-1,1)  &
-         !           + tmask(ji-1,jj+1,1) + tmask(ji+1,jj+1,1) ) 
-         !   ENDDO
-         !ENDDO
-
      ENDIF
-     !PRINT*,"min max zmask ",MINVAL(zmask),MAXVAL(zmask)
 
      !debug
@@ -148,6 +127,6 @@
 
      !loop on the mesh 
-     DO jj=2,nlcj
-        DO ji=2,nlci
+     DO jj=2,jpj
+        DO ji=2,jpi
            !-----------------------------------------------------------
            !For each cell of the mesh, we find the intersection between 
@@ -176,66 +155,66 @@
            zdeltaj= gphif(ji,jj-1) - gphif(ji-1,jj-1) !zdeltaj=[AD]
            
-           IF (ABS(zdeltai) .LE. 2.*zdistmesh .OR.  & 
-               ABS(zdeltaj) .LE. 2.*zdistmesh) THEN           
+           IF( ABS(zdeltai) .LE. 2.*zdistmesh .OR.  & 
+               ABS(zdeltaj) .LE. 2.*zdistmesh )THEN           
   
-               !intersection section/[AB]=?
-               coordTemp     = intersec(sec,coord_a,coord_b)     !compute intersection 
-               coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
-               coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
- 
-               IF(coordTemp%lon .NE. -9999) THEN
-                   IF(nb_inmesh+1 .GT. nb_point_max) THEN
-                       PRINT*,"WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh
-                   ELSE    
-                       nb_inmesh=nb_inmesh+1              
-                       coordSec(nb_inmesh) = coordTemp    !store the intersection's coordinates
-
-                       !We need to know if the section crosses the overlapping band.
-
-                       !Fist we look if there is an intersection mesh/section 
-                       !just on the left of the overlapping band: 
-                       IF(coordTemp%lon .GT. glamf(1,1)-5  .AND. coordTemp%lon .LT. glamf(1,1)) & 
+              !intersection section/[AB]=?
+              coordTemp     = intersec(sec,coord_a,coord_b)     !compute intersection 
+              coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
+              coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
+ 
+              IF( coordTemp%lon .NE. -9999 )THEN
+                 IF( nb_inmesh+1 .GT. nb_point_max )THEN
+                    PRINT*,"WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh
+                 ELSE    
+                    nb_inmesh=nb_inmesh+1              
+                    coordSec(nb_inmesh) = coordTemp    !store the intersection's coordinates
+
+                    !We need to know if the section crosses the overlapping band.
+
+                    !Fist we look if there is an intersection mesh/section 
+                    !just on the left of the overlapping band: 
+                    IF( coordTemp%lon .GT. glamf(1,1)-5  .AND. coordTemp%lon .LT. glamf(1,1) ) & 
+                       & ll_overlap_sec_left  = .TRUE.
+                    !And we look if there is an intersection mesh/section 
+                    !just on the right of the overlapping band:
+                    IF( coordTemp%lon .GT. glamf(jpi,1) .AND. coordTemp%lon .LT. glamf(1,1)+5 ) &
+                       & ll_overlap_sec_right = .TRUE.
+                 ENDIF
+              ENDIF
+              
+              !intersection section/[AD]=?
+              coordTemp=intersec(sec,coord_a,coord_d)       !compute intersection
+              coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
+              coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
+
+              IF( coordTemp%lon .NE. -9999 )THEN
+                 IF( nb_inmesh+1 .GT. nb_point_max )THEN
+                    PRINT*, "WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh
+                 ELSE  
+                    nb_inmesh=nb_inmesh+1
+                    coordSec(nb_inmesh)=coordTemp
+                   
+                    !We need to know if the section crosses the overlapping band:
+                    !same test as above 
+                    IF( coordTemp%lon .GE. glamf(1,1)-3  .AND. coordTemp%lon .LE. glamf(1,1) ) &
                          & ll_overlap_sec_left  = .TRUE.
-                       !And we look if there is an intersection mesh/section 
-                       !just on the right of the overlapping band:
-                       IF(coordTemp%lon .GT. glamf(nlci,1) .AND. coordTemp%lon .LT. glamf(1,1)+5) &
-                         & ll_overlap_sec_right = .TRUE.
-                   ENDIF
-               ENDIF
-              
-               !intersection section/[AD]=?
-               coordTemp=intersec(sec,coord_a,coord_d)       !compute intersection
-               coordTemp%lon = coordTemp%lon !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
-               coordTemp%lat = coordTemp%lat !* zmask(ji,jj)-9999.*(1-zmask(ji,jj))
-
-               IF(coordTemp%lon .NE. -9999) THEN
-                   IF(nb_inmesh+1 .GT. nb_point_max) THEN
-                       PRINT*, "WARNING diadct: nb_point_max needs to be greater than ", nb_inmesh
-                   ELSE  
-                       nb_inmesh=nb_inmesh+1
-                       coordSec(nb_inmesh)=coordTemp
-                   
-                       !We need to know if the section crosses the overlapping band:
-                       !same test as above 
-                       IF(coordTemp%lon .GE. glamf(1,1)-3  .AND. coordTemp%lon .LE. glamf(1,1)) &
-                         & ll_overlap_sec_left  = .TRUE.
-                       IF(coordTemp%lon .GE. glamf(nlci,1) .AND. coordTemp%lon .LE. glamf(nlci,1)+3) & 
+                    IF( coordTemp%lon .GE. glamf(jpi,1) .AND. coordTemp%lon .LE. glamf(jpi,1)+3) & 
                          & ll_overlap_sec_right = .TRUE. 
-                   ENDIF
-               ENDIF
- 
-          ENDIF
+                 ENDIF
+              ENDIF
+ 
+           ENDIF
          
-          !We need to know if the domain crosses the date line:
-          !Fist, we search a mesh point that is just one the left of date line:
-          IF( glamf(ji-1,jj-1) .GT.  175 .AND. glamf(ji-1,jj-1) .LT.  180 ) &
-            & ll_date_domain_left = .TRUE.
-          !And we search a mesh point that is just one the right of date line:
-          IF( glamf(ji-1,jj-1) .GT. -180 .AND. glamf(ji-1,jj-1) .LT. -175 ) &  
-            & ll_date_domain_right = .TRUE.
- 
-         !End of the loop on the mesh
-       ENDDO
-     ENDDO
+           !We need to know if the domain crosses the date line:
+           !Fist, we search a mesh point that is just one the left of date line:
+           IF( glamf(ji-1,jj-1) .GT.  175 .AND. glamf(ji-1,jj-1) .LT.  180 ) &
+              & ll_date_domain_left = .TRUE.
+           !And we search a mesh point that is just one the right of date line:
+           IF( glamf(ji-1,jj-1) .GT. -180 .AND. glamf(ji-1,jj-1) .LT. -175 ) &  
+              & ll_date_domain_right = .TRUE.
+ 
+        ENDDO
+     ENDDO !End of the loop on the mesh
+
  
      !Crossing section/overlapping band (we need to know it for later):
@@ -266,22 +245,26 @@
 
      IF( sec%ll_date_line .AND. ll_date_domain )THEN
+
         !we add 360° to negative longitudes to have a good classification
-         DO jseg=1,nb_inmesh
-              IF( coordSec(jseg)%lon .LT. 0 ) coordSec(jseg)%lon=coordSec(jseg)%lon+360.
-         ENDDO
-         IF ( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon ) THEN
-              CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh)
-         ELSE
-              CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh)
-         ENDIF
-         DO jseg=1,nb_inmesh
-              IF( coordSec(jseg)%lon .GT. 180 ) coordSec(jseg)%lon=coordSec(jseg)%lon-360.
-         ENDDO
+        DO jseg=1,nb_inmesh
+           IF( coordSec(jseg)%lon .LT. 0 ) coordSec(jseg)%lon=coordSec(jseg)%lon+360.
+        ENDDO
+        IF( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon ) THEN
+           CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh)
+        ELSE
+           CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh)
+        ENDIF
+        DO jseg=1,nb_inmesh
+           IF( coordSec(jseg)%lon .GT. 180 ) coordSec(jseg)%lon=coordSec(jseg)%lon-360.
+        ENDDO
+
      ELSE     
-         IF ( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon ) THEN
-              CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh)
-         ELSE
-              CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh)
-         ENDIF
+
+        IF( sec%coordSec(1)%lon .NE. sec%coordSec(2)%lon )THEN
+           CALL qcksrt(coordSec(:)%lon,coordSec(:)%lat,nb_inmesh)
+        ELSE
+           CALL qcksrt(coordSec(:)%lat,coordSec(:)%lon,nb_inmesh)
+        ENDIF
+
      ENDIF
 
@@ -300,49 +283,49 @@
 
      IF( nb_inmesh .ne. 0 )THEN
-         coordFirst       = coordSec(1)
-         coordLast        = coordSec(nb_inmesh) 
-         sec%nb_point     = nb_inmesh
-         sec%listPoint(1) = POINT_SECTION(-1,-1)
-         zdistante        = 1000.
-         zdistante2       = 1000.
-
-         !First, we find the point of the mesh that is the closest 
-         !to the first intersection section/mesh (=coordFirst=coordSec(1)):
-         !this point will be called sec%listPoint(1).
-         !Then, we find the point of the mesh that is the closest
-         !to the last intersection section/mesh (coordLast=coordSec(nb_inmesh))
-         !this point will be called endingPoint.
-
-         DO jj=1,nlej
-            DO ji=1,nlei  
-               coord_t=COORD_SECTION(glamf(ji,jj),gphif(ji,jj))
-               zdistFirst = distance(coord_t,coordFirst)
-               zdistLast = distance(coord_t,coordLast)
-               IF( zdistFirst .LT. zdistmesh .AND. zdistFirst .LT. zdistante )THEN
-                    sec%listPoint(1) = POINT_SECTION(ji,jj)
-                    zdistante=zdistFirst
-               ENDIF
-               IF( zdistLast .LT. zdistmesh .AND. zdistLast .LT. zdistante2 )THEN
-                    endingPoint = POINT_SECTION(ji,jj)
-                    zdistante2=zdistLast
-               ENDIF
-            ENDDO
-         ENDDO
-
-         IF( sec%listPoint(1)%I == endingPoint%I .AND. sec%listPoint(1)%J == endingPoint%J )THEN
-            sec%listPoint(1) = POINT_SECTION(-1,-1)
-            endingPoint      = POINT_SECTION(-1,-1)
-            coordFirst       = coordSec(1)
-            coordLast        = coordSec(2)
-            sec%nb_point     = 0
-         ENDIF
+        coordFirst       = coordSec(1)
+        coordLast        = coordSec(nb_inmesh) 
+        sec%nb_point     = nb_inmesh
+        sec%listPoint(1) = POINT_SECTION(-1,-1)
+        zdistante        = 1000.
+        zdistante2       = 1000.
+
+        !First, we find the point of the mesh that is the closest 
+        !to the first intersection section/mesh (=coordFirst=coordSec(1)):
+        !this point will be called sec%listPoint(1).
+        !Then, we find the point of the mesh that is the closest
+        !to the last intersection section/mesh (coordLast=coordSec(nb_inmesh))
+        !this point will be called endingPoint.
+
+        DO jj=1,jpj
+           DO ji=1,jpi  
+              coord_t=COORD_SECTION(glamf(ji,jj),gphif(ji,jj))
+              zdistFirst = distance(coord_t,coordFirst)
+              zdistLast = distance(coord_t,coordLast)
+              IF( zdistFirst .LT. zdistmesh .AND. zdistFirst .LT. zdistante )THEN
+                 sec%listPoint(1) = POINT_SECTION(ji,jj)
+                 zdistante=zdistFirst
+              ENDIF
+              IF( zdistLast .LT. zdistmesh .AND. zdistLast .LT. zdistante2 )THEN
+                 endingPoint = POINT_SECTION(ji,jj)
+                 zdistante2=zdistLast
+              ENDIF
+           ENDDO
+        ENDDO
+
+        IF( sec%listPoint(1)%I == endingPoint%I .AND. sec%listPoint(1)%J == endingPoint%J )THEN
+           sec%listPoint(1) = POINT_SECTION(-1,-1)
+           endingPoint      = POINT_SECTION(-1,-1)
+           coordFirst       = coordSec(1)
+           coordLast        = coordSec(2)
+           sec%nb_point     = 0
+        ENDIF
 
      ELSE
-         !If there is no intersection section/mesh
-         sec%listPoint(1) = POINT_SECTION(-1,-1)
-         endingPoint      = POINT_SECTION(-1,-1)
-         coordFirst       = coordSec(1)
-         coordLast        = coordSec(2)
-         sec%nb_point     = 0
+        !If there is no intersection section/mesh
+        sec%listPoint(1) = POINT_SECTION(-1,-1)
+        endingPoint      = POINT_SECTION(-1,-1)
+        coordFirst       = coordSec(1)
+        coordLast        = coordSec(2)
+        sec%nb_point     = 0
      ENDIF 
 
@@ -380,294 +363,295 @@
      IF( nb_inmesh .NE. 0 )THEN
 
-         !The serie of mesh's points that form the section will 'link' 
-         !sec%listPoint(1) to endingPoint: it will be stored in 
-         !sec%listPoint(jseg)
-         !
-         !We take place on the fist point (sec%listPoint(1)) 
-         ! a.  We find the 4 adjacent points (North, South, East, West)
-         ! b.  Compute distance between current point and endingPoint
-         ! c.  Compute distance between the 4 adjacent points and endingPoint
-         ! d.  Select the points which are closer to end-point than current point
-         ! e.1 If at least one point is selected, select the point which is closest to original section among selected points 
-         ! e.2 If no point is selected, select the point which is the closest to end-point 
-         ! f. save next point and direction of velocity.
-         ! g. Save nextPoint and go to nextPoint
-         !
-         !We get out of this loop if:
-         !    - we are on endingPoint
-         !    - the number of points (jseg) that link sec%listPoint(1) to endingPoint is 
-         !      twice greater than number of section/mesh intersection (nb_inmesh):
-         !      it could be possible if thr algorithm can't link endingPoint (bug).
-
-         !initialize distnew value (with distance between section's extremities) 
-         zdistnew  = distance(coordFirst,coordLast,sec%ll_date_line) 
-         prevPoint  = POINT_SECTION(0,0)
-         jseg       = 1
-
-         DO WHILE ( (  sec%listPoint(jseg)%I .NE.  endingPoint%I    &
-                  .OR. sec%listPoint(jseg)%J .NE. endingPoint%J   ) &
-                  .AND. jseg .LT. 500 .AND. sec%listPoint(jseg)%I .GT. 0  )         
+        !The serie of mesh's points that form the section will 'link' 
+        !sec%listPoint(1) to endingPoint: it will be stored in 
+        !sec%listPoint(jseg)
+        !
+        !We take place on the fist point (sec%listPoint(1)) 
+        ! a.  We find the 4 adjacent points (North, South, East, West)
+        ! b.  Compute distance between current point and endingPoint
+        ! c.  Compute distance between the 4 adjacent points and endingPoint
+        ! d.  Select the points which are closer to end-point than current point
+        ! e.1 If at least one point is selected, select the point which is closest to original section among selected points 
+        ! e.2 If no point is selected, select the point which is the closest to end-point 
+        ! f. save next point and direction of velocity.
+        ! g. Save nextPoint and go to nextPoint
+        !
+        !We get out of this loop if:
+        !    - we are on endingPoint
+        !    - the number of points (jseg) that link sec%listPoint(1) to endingPoint is 
+        !      twice greater than number of section/mesh intersection (nb_inmesh):
+        !      it could be possible if thr algorithm can't link endingPoint (bug).
+
+        !initialize distnew value (with distance between section's extremities) 
+        zdistnew  = distance(coordFirst,coordLast,sec%ll_date_line) 
+        prevPoint  = POINT_SECTION(0,0)
+        jseg       = 1
+
+        DO WHILE ( (  sec%listPoint(jseg)%I .NE.  endingPoint%I    &
+                 .OR. sec%listPoint(jseg)%J .NE. endingPoint%J   ) &
+                 .AND. jseg .LT. 500 .AND. sec%listPoint(jseg)%I .GT. 0  )         
    
-            ! a. find the 4 adjacent points (North, South, East, West)
-            !---------------------------------------------------------
-            SouthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J-1)
-            NorthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
-            WestPoint  = POINT_SECTION(sec%listPoint(jseg)%I-1,sec%listPoint(jseg)%J)
-            EstPoint   = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
-
-            !debug
-            CALL write_debug(jsec,"---------------")
-            WRITE(cltmp,100)'Current points: ',sec%listPoint(jseg), & 
-                          glamf(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J), &
-                          gphif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)
-            CALL write_debug(jsec,cltmp)
-            CALL write_debug(jsec,"E/W/N/S points")
-            WRITE(cltmp,101)glamf(  EstPoint%I,EstPoint%J)  ,gphif(  EstPoint%I,  EstPoint%J), &
-                            glamf( WestPoint%I,WestPoint%J) ,gphif( WestPoint%I, WestPoint%J), &
-                            glamf(NorthPoint%I,NorthPoint%J),gphif(NorthPoint%I,NorthPoint%J) ,&
-                            glamf(SouthPoint%I,SouthPoint%J),gphif(SouthPoint%I,SouthPoint%J)
-            CALL write_debug(jsec,cltmp)
-            WRITE(cltmp,102)EstPoint,WestPoint,NorthPoint,SouthPoint
-            CALL write_debug(jsec,cltmp)
+           ! a. find the 4 adjacent points (North, South, East, West)
+           !---------------------------------------------------------
+           SouthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J-1)
+           NorthPoint = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
+           WestPoint  = POINT_SECTION(sec%listPoint(jseg)%I-1,sec%listPoint(jseg)%J)
+           EstPoint   = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
+
+           !debug
+           CALL write_debug(jsec,"---------------")
+           WRITE(cltmp,100)'Current points: ',sec%listPoint(jseg), & 
+                         glamf(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J), &
+                         gphif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)
+           CALL write_debug(jsec,cltmp)
+           CALL write_debug(jsec,"E/W/N/S points")
+           WRITE(cltmp,101)glamf(  EstPoint%I,EstPoint%J)  ,gphif(  EstPoint%I,  EstPoint%J), &
+                           glamf( WestPoint%I,WestPoint%J) ,gphif( WestPoint%I, WestPoint%J), &
+                           glamf(NorthPoint%I,NorthPoint%J),gphif(NorthPoint%I,NorthPoint%J) ,&
+                           glamf(SouthPoint%I,SouthPoint%J),gphif(SouthPoint%I,SouthPoint%J)
+           CALL write_debug(jsec,cltmp)
+           WRITE(cltmp,102)EstPoint,WestPoint,NorthPoint,SouthPoint
+           CALL write_debug(jsec,cltmp)
 
 100 FORMAT ( A15,2(i4.4," "),2(f7.3," ") )
 101 FORMAT ( "E ",2(f7.3," "),"W ",2(f7.3," "),"N ",2(f7.3," "),"S ",2(f7.3," "))
-102 FORMAT ( "E ",i4.4,' 'i4.4,"//W ",i4.4,' ',i4.4,"//N ",i4.4,' ',i4.4,"//S ",i4.4,' ',i4.4 )
+102 FORMAT ( "E ",i4.4,' ',i4.4,"//W ",i4.4,' ',i4.4,"//N ",i4.4,' ',i4.4,"//S ",i4.4,' ',i4.4 )
 
                
-            !we  are on end-point
-            !--------------------
-            IF(      SouthPoint%I==endingPoint%I .AND. SouthPoint%J==endingPoint%J )THEN 
-                jseg = jseg+1 ; sec%listPoint(jseg) = SouthPoint
-            ELSE IF( NorthPoint%I==endingPoint%I .AND. NorthPoint%J==endingPoint%J )THEN
-                jseg = jseg+1 ; sec%listPoint(jseg) = NorthPoint
-            ELSE IF(  WestPoint%I==endingPoint%I .AND.  WestPoint%J==endingPoint%J )THEN
-                jseg = jseg+1 ; sec%listPoint(jseg) = WestPoint
-            ELSE IF(   EstPoint%I==endingPoint%I .AND.   EstPoint%J==endingPoint%J )THEN
-                jseg = jseg+1 ; sec%listPoint(jseg) = EstPoint
-
-            ELSE
-            !we  are NOT on end-point
-            !------------------------
-
-               ! b. distance between current point and endingPoint
-               !--------------------------------------------------
-               zdistante = zdistnew
-
-               ! c. compute distance between the 4 adjacent points and endingPoint
-               !------------------------------------------------------------------
-               ! BE CARREFUL! When the domain crosses the date line (ll_date_domain):
-               ! When we will compute distances between W/E/S/N points and endingPoint,
-               ! we have to check if theses 4 lines crosses the date line
-               ! (test: ABS(coordTemp%lon - coordLast%lon).GT. 180)
-               ! If it's true,we have to add 360° to coordLast%long to compute the 
-               ! good distance through the date line and not through the complementary
-               ! in the mesh.
+           !we  are on end-point
+           !--------------------
+           IF(      SouthPoint%I==endingPoint%I .AND. SouthPoint%J==endingPoint%J )THEN 
+               jseg = jseg+1 ; sec%listPoint(jseg) = SouthPoint
+           ELSE IF( NorthPoint%I==endingPoint%I .AND. NorthPoint%J==endingPoint%J )THEN
+               jseg = jseg+1 ; sec%listPoint(jseg) = NorthPoint
+           ELSE IF(  WestPoint%I==endingPoint%I .AND.  WestPoint%J==endingPoint%J )THEN
+               jseg = jseg+1 ; sec%listPoint(jseg) = WestPoint
+           ELSE IF(   EstPoint%I==endingPoint%I .AND.   EstPoint%J==endingPoint%J )THEN
+               jseg = jseg+1 ; sec%listPoint(jseg) = EstPoint
+
+           ELSE
+           !we  are NOT on end-point
+           !------------------------
+
+              ! b. distance between current point and endingPoint
+              !--------------------------------------------------
+              zdistante = zdistnew
+
+              ! c. compute distance between the 4 adjacent points and endingPoint
+              !------------------------------------------------------------------
+              ! BE CARREFUL! When the domain crosses the date line (ll_date_domain):
+              ! When we will compute distances between W/E/S/N points and endingPoint,
+              ! we have to check if theses 4 lines crosses the date line
+              ! (test: ABS(coordTemp%lon - coordLast%lon).GT. 180)
+              ! If it's true,we have to add 360° to coordLast%long to compute the 
+              ! good distance through the date line and not through the complementary
+              ! in the mesh.
        
-               ! c.1 compute distWest: distance between west point and endingPoint
-               !----------------------
-               zdistWest2 = 99999999.9 ;  zdistWest3 = 99999999.9 
-               IF( sec%listPoint(jseg)%I .NE. 1 )THEN
-                  !When we are on the west side of the mesh (i=1),we can't go to the west.
-                  coordTemp = pointToCoordF(WestPoint) 
-                  ll_test = .FALSE.
-                  IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test = .TRUE.
-                  zdistWest2  = distance2(pointToCoordF(WestPoint) ,coordLast ,ll_test)
-               ENDIF
-
-               ! c.2 compute distEst: distance between east point and endingPoint
-               !---------------------
-               zdistEst2 = 99999999.9 ;  zdistEst3 = 99999999.9
-               IF( sec%listPoint(jseg)%I .EQ. nlci )THEN
-                  !We test if the current point is on the east side of the mesh
-                  ! The method is done such as we go toward east to link
-                  ! sec%listPoint(1) to endingPoint.
-                  ! So, if the section crosses the overlapping band (ll_overlap_sec=T),
-                  ! we won't have to stop if the current point is on the EAST side of the mesh:
-                  ! we have to follow the construction of the section on the
-                  ! WEST side of the mesh
-                  IF( ll_overlap_sec  )THEN
-                     !section crosses the overlapping band 
-                     !So EstPoint is on the west side of the mesh
-                     EstPoint = POINT_SECTION(3,sec%listPoint(jseg)%J)
-                     zdistEst2= distance2(pointToCoordF(EstPoint)  ,coordLast ,.FALSE.)
-                  ENDIF 
-               ELSE
-                  coordTemp = pointToCoordF(EstPoint) 
-                  ll_test = .FALSE.
-                  IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
-                  zdistEst2 = distance2(pointToCoordF(EstPoint)  ,coordLast ,ll_test )
-               ENDIF
-
-               ! c.3 compute distSouth: distance between south point and endingPoint
-               !-----------------------
-               zdistSouth2 = 99999999.9 ; zdistSouth3 = 99999999.9
-               IF( sec%listPoint(jseg)%J .NE. 1 )THEN
-                  !When we are on the south side of the mesh (j=1),we can't go to the south. 
-                  coordTemp=pointToCoordF(SouthPoint)
-                  ll_test = .FALSE.
-                  IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
-                     zdistSouth2 = distance2(pointToCoordF(SouthPoint),coordlast ,ll_test )
-                  ENDIF
-
-               ! c.4 compute distNorth: distance between north and endingPoint
-               !-----------------------
-               zdistNorth2 = 99999999.9 ; zdistNorth3 = 99999999.9 
-               IF( sec%listPoint(jseg)%J .NE. nlcj )THEN
-                  !When we are on the north side of the mesh (j=nlcj),we can't go to the south. 
-                  coordTemp=pointToCoordF(NorthPoint)
-                  ll_test = .FALSE.
-                  IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
-                  zdistNorth2 = distance2(pointToCoordF(NorthPoint),coordlast ,ll_test )
-               ENDIF
-
-               ! d. select the points which are closer to end-point than current point
-               !----------------------------------------------------------------------
-               zdistref=distance2(pointToCoordF(sec%listPoint(jseg)),coordlast ,ll_test )
-               WRITE(cltmp,'( A56,f10.3 )' )'distance between actual point and last point: zdistref = ',zdistref
-               CALL write_debug(jsec,cltmp)
-               lest   = .FALSE. ; IF( zdistEst2   .LE. zdistref ) lest  = .TRUE.
-               lwest  = .FALSE. ; IF( zdistwest2  .LE. zdistref ) lwest = .TRUE.
-               lnorth = .FALSE. ; IF( zdistnorth2 .LE. zdistref ) lnorth= .TRUE.
-               lsouth = .FALSE. ; IF( zdistsouth2 .LE. zdistref ) lsouth= .TRUE.
-
-               !debug
-               IF( .NOT. lest   )CALL write_debug(jsec,'Est   point eliminated')
-               IF( .NOT. lwest  )CALL write_debug(jsec,'West  point eliminated')
-               IF( .NOT. lnorth )CALL write_debug(jsec,'North point eliminated')
-               IF( .NOT. lsouth )CALL write_debug(jsec,'South point eliminated')
-
-               l_oldmethod = .FALSE.
-
-               IF( ( COUNT((/lest/))+COUNT((/lwest/))+COUNT((/lnorth/))+COUNT((/lsouth/)) ) .NE. 0 )THEN
-
-                  ! e.1 If at least one point is selected, select the point 
-                  !     which is the closest to original section among selected points 
-                  !-------------------------------------------------------------------
-
-                  zdistWest3  = 9999999.9
-                  IF( lwest )zdistWest3  = &
-                     distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(WestPoint) ,pointToCoordF(endingPoint) ,lkdebug )
-                  zdistEst3   = 9999999.9
-                  IF( lest )zdistEst3   =  &
-                     distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(EstPoint)  ,pointToCoordF(endingPoint) ,lkdebug )
-                  zdistSouth3 = 9999999.9
-                  IF( lsouth )zdistSouth3 = &
-                     distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(SouthPoint),pointToCoordF(endingPoint) ,lkdebug )
-                  zdistNorth3 = 9999999.9
-                  IF( lnorth )zdistNorth3 = &
-                     distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(NorthPoint),pointToCoordF(endingPoint) ,lkdebug )
-
-                  zdistEst3   = zdistEst3   + (1-COUNT((/lest/))  )*9999999.9
-                  zdistWest3  = zdistWest3  + (1-COUNT((/lwest/)) )*9999999.9
-                  zdistNorth3 = zdistNorth3 + (1-COUNT((/lnorth/)))*9999999.9
-                  zdistSouth3 = zdistSouth3 + (1-COUNT((/lsouth/)))*9999999.9
-
-                  zdistnew = MIN(zdistEst3,zdistWest3,zdistnorth3,zdistSouth3)
-
-               ELSE 
-
-                  ! e.2 If no point is selected, select the point which is the closest to end-point 
-                  !--------------------------------------------------------------------------------
-                  l_oldmethod = .TRUE.
-
-                  !debug
-                  WRITE(cltmp,'(A30,i3.3)' )'SEARCH NEW POINT WITH OLD METHOD: ',jsec
-                  CALL write_debug(jsec,cltmp)
+              ! c.1 compute distWest: distance between west point and endingPoint
+              !----------------------
+              zdistWest2 = 99999999.9 ;  zdistWest3 = 99999999.9 
+              IF( sec%listPoint(jseg)%I .NE. 1 )THEN
+                 !When we are on the west side of the mesh (i=1),we can't go to the west.
+                 coordTemp = pointToCoordF(WestPoint) 
+                 ll_test = .FALSE.
+                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test = .TRUE.
+                 zdistWest2  = distance2(pointToCoordF(WestPoint) ,coordLast ,ll_test)
+              ENDIF
+
+              ! c.2 compute distEst: distance between east point and endingPoint
+              !---------------------
+              zdistEst2 = 99999999.9 ;  zdistEst3 = 99999999.9
+              IF( sec%listPoint(jseg)%I .EQ. jpi )THEN
+                 !We test if the current point is on the east side of the mesh
+                 ! The method is done such as we go toward east to link
+                 ! sec%listPoint(1) to endingPoint.
+                 ! So, if the section crosses the overlapping band (ll_overlap_sec=T),
+                 ! we won't have to stop if the current point is on the EAST side of the mesh:
+                 ! we have to follow the construction of the section on the
+                 ! WEST side of the mesh
+                 IF( ll_overlap_sec  )THEN
+                    !section crosses the overlapping band 
+                    !So EstPoint is on the west side of the mesh
+                    EstPoint = POINT_SECTION(3,sec%listPoint(jseg)%J)
+                    zdistEst2= distance2(pointToCoordF(EstPoint)  ,coordLast ,.FALSE.)
+                 ENDIF 
+              ELSE
+                 coordTemp = pointToCoordF(EstPoint) 
+                 ll_test = .FALSE.
+                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
+                 zdistEst2 = distance2(pointToCoordF(EstPoint)  ,coordLast ,ll_test )
+              ENDIF
+
+              ! c.3 compute distSouth: distance between south point and endingPoint
+              !-----------------------
+              zdistSouth2 = 99999999.9 ; zdistSouth3 = 99999999.9
+              IF( sec%listPoint(jseg)%J .NE. 1 )THEN
+                 !When we are on the south side of the mesh (j=1),we can't go to the south. 
+                 coordTemp=pointToCoordF(SouthPoint)
+                 ll_test = .FALSE.
+                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
+                    zdistSouth2 = distance2(pointToCoordF(SouthPoint),coordlast ,ll_test )
+              ENDIF
+
+              ! c.4 compute distNorth: distance between north and endingPoint
+              !-----------------------
+              zdistNorth2 = 99999999.9 ; zdistNorth3 = 99999999.9 
+              IF( sec%listPoint(jseg)%J .NE. jpj )THEN
+                 !When we are on the north side of the mesh (j=jpj),we can't go to the south. 
+                 coordTemp=pointToCoordF(NorthPoint)
+                 ll_test = .FALSE.
+                 IF( ll_date_domain .AND. ABS(coordTemp%lon - coordLast%lon).GT. 180 )ll_test= .TRUE.
+                 zdistNorth2 = distance2(pointToCoordF(NorthPoint),coordlast ,ll_test )
+              ENDIF
+
+              ! d. select the points which are closer to end-point than current point
+              !----------------------------------------------------------------------
+              zdistref=distance2(pointToCoordF(sec%listPoint(jseg)),coordlast ,ll_test )
+              WRITE(cltmp,'( A56,f10.3 )' )'distance between actual point and last point: zdistref = ',zdistref
+              CALL write_debug(jsec,cltmp)
+              lest   = .FALSE. ; IF( zdistEst2   .LE. zdistref ) lest  = .TRUE.
+              lwest  = .FALSE. ; IF( zdistwest2  .LE. zdistref ) lwest = .TRUE.
+              lnorth = .FALSE. ; IF( zdistnorth2 .LE. zdistref ) lnorth= .TRUE.
+              lsouth = .FALSE. ; IF( zdistsouth2 .LE. zdistref ) lsouth= .TRUE.
+
+              !debug
+              IF( .NOT. lest   )CALL write_debug(jsec,'Est   point eliminated')
+              IF( .NOT. lwest  )CALL write_debug(jsec,'West  point eliminated')
+              IF( .NOT. lnorth )CALL write_debug(jsec,'North point eliminated')
+              IF( .NOT. lsouth )CALL write_debug(jsec,'South point eliminated')
+
+              l_oldmethod = .FALSE.
+
+              IF( ( COUNT((/lest/))+COUNT((/lwest/))+COUNT((/lnorth/))+COUNT((/lsouth/)) ) .NE. 0 )THEN
+
+                 ! e.1 If at least one point is selected, select the point 
+                 !     which is the closest to original section among selected points 
+                 !-------------------------------------------------------------------
+
+                 zdistWest3  = 9999999.9
+                 IF( lwest )zdistWest3  = &
+                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(WestPoint) ,pointToCoordF(endingPoint) ,lkdebug )
+                 zdistEst3   = 9999999.9
+                 IF( lest )zdistEst3   =  &
+                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(EstPoint)  ,pointToCoordF(endingPoint) ,lkdebug )
+                 zdistSouth3 = 9999999.9
+                 IF( lsouth )zdistSouth3 = &
+                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(SouthPoint),pointToCoordF(endingPoint) ,lkdebug )
+                 zdistNorth3 = 9999999.9
+                 IF( lnorth )zdistNorth3 = &
+                    distance3(pointToCoordF(sec%listPoint(1)),pointToCoordF(NorthPoint),pointToCoordF(endingPoint) ,lkdebug )
+
+                 zdistEst3   = zdistEst3   + (1-COUNT((/lest/))  )*9999999.9
+                 zdistWest3  = zdistWest3  + (1-COUNT((/lwest/)) )*9999999.9
+                 zdistNorth3 = zdistNorth3 + (1-COUNT((/lnorth/)))*9999999.9
+                 zdistSouth3 = zdistSouth3 + (1-COUNT((/lsouth/)))*9999999.9
+
+                 zdistnew = MIN(zdistEst3,zdistWest3,zdistnorth3,zdistSouth3)
+
+              ELSE 
+
+                 ! e.2 If no point is selected, select the point which is the closest to end-point 
+                 !--------------------------------------------------------------------------------
+                 l_oldmethod = .TRUE.
+
+                 !debug
+                 WRITE(cltmp,'(A30,i3.3)' )'SEARCH NEW POINT WITH OLD METHOD: ',jsec
+                 CALL write_debug(jsec,cltmp)
                    
-                  ! on passe à l'ancienne methode: le point parmies les 4 pts (NSWE)  qui se rapproche
-                  ! le + du dernier pt
-                  !----------------------------- 
-                  !be carreful! we can't go backward.
-
-                  zdistNorth = zdistNorth2    ; zdistSouth = zdistSouth2
-                  zdistEst   = zdistEst2      ; zdistWest  = zdistWest2 
-
-                  IF(     prevPoint%I .EQ. NorthPoint%I .AND. prevPoint%J .EQ. NorthPoint%J) THEN
-                      zdistnew = MIN(zdistEst,zdistWest,zdistSouth)
-                  ELSE IF(prevPoint%I .EQ. SouthPoint%I .AND. prevPoint%J .EQ. SouthPoint%J) THEN
-                     zdistnew = MIN(zdistEst,zdistWest,zdistNorth)
-                  ELSE IF(prevPoint%I .EQ. WestPoint%I  .AND. prevPoint%J .EQ. WestPoint%J ) THEN
-                     zdistnew = MIN(zdistEst,zdistNorth,zdistSouth)
-                  ELSE IF(prevPoint%I .EQ. EstPoint%I   .AND. prevPoint%J .EQ. EstPoint%J  ) THEN
-                     zdistnew = MIN(zdistWest,zdistNorth,zdistSouth)
-                  ELSE 
-                     zdistnew = MIN(zdistEst,zdistWest,zdistNorth,zdistSouth)
-                  ENDIF              
-
-               ENDIF
-
-               !debug
-               WRITE(cltmp,'(A11, f8.3)')'zdistref = ',zdistref
-               CALL write_debug(jsec,cltmp)
-               WRITE(cltmp, 103         )'distance2 :',zdistEst2,zdistWest2,zdistNorth2,zdistSouth2 
-               CALL write_debug(jsec,cltmp)
-               WRITE(cltmp, 103         )'distance3 :',zdistEst3,zdistWest3,zdistNorth3,zdistSouth3
-               CALL write_debug(jsec,cltmp)
-               WRITE(cltmp,'(A11, f8.3)')"zdistnew = ",zdistnew
-               CALL write_debug(jsec,cltmp)
+                 ! on passe à l'ancienne methode: le point parmies les 4 pts (NSWE)  qui se rapproche
+                 ! le + du dernier pt
+                 !----------------------------- 
+                 !be carreful! we can't go backward.
+
+                 zdistNorth = zdistNorth2    ; zdistSouth = zdistSouth2
+                 zdistEst   = zdistEst2      ; zdistWest  = zdistWest2 
+
+                 IF(     prevPoint%I .EQ. NorthPoint%I .AND. prevPoint%J .EQ. NorthPoint%J) THEN
+                     zdistnew = MIN(zdistEst,zdistWest,zdistSouth)
+                 ELSE IF(prevPoint%I .EQ. SouthPoint%I .AND. prevPoint%J .EQ. SouthPoint%J) THEN
+                    zdistnew = MIN(zdistEst,zdistWest,zdistNorth)
+                 ELSE IF(prevPoint%I .EQ. WestPoint%I  .AND. prevPoint%J .EQ. WestPoint%J ) THEN
+                    zdistnew = MIN(zdistEst,zdistNorth,zdistSouth)
+                 ELSE IF(prevPoint%I .EQ. EstPoint%I   .AND. prevPoint%J .EQ. EstPoint%J  ) THEN
+                    zdistnew = MIN(zdistWest,zdistNorth,zdistSouth)
+                 ELSE 
+                    zdistnew = MIN(zdistEst,zdistWest,zdistNorth,zdistSouth)
+                 ENDIF              
+
+              ENDIF
+
+              !debug
+              WRITE(cltmp,'(A11, f8.3)')'zdistref = ',zdistref
+              CALL write_debug(jsec,cltmp)
+              WRITE(cltmp, 103         )'distance2 :',zdistEst2,zdistWest2,zdistNorth2,zdistSouth2 
+              CALL write_debug(jsec,cltmp)
+              WRITE(cltmp, 103         )'distance3 :',zdistEst3,zdistWest3,zdistNorth3,zdistSouth3
+              CALL write_debug(jsec,cltmp)
+              WRITE(cltmp,'(A11, f8.3)')"zdistnew = ",zdistnew
+              CALL write_debug(jsec,cltmp)
 
 103 FORMAT (A12,"E",f12.3," W",f12.3," N",f12.3," S",f12.3)
 
-               !f. save next point and direction of velocity.
-               !--------------------------------------------- 
-               !nextPoint will be the one which is the closest to endingPoint.
-               !sec%direction will be direction between current point and nextPoint:
-               !It will be used to compute velocity through the segment [CurrentPoint,nextPoint}:
-               !                 
-               !A:Current Point    NorthPoint(I,J+1)   Nextpoint=NorthPoint(I,J+1) => sec%direction=3
-               !                   |                   Nextpoint=SouthPoint(I,J-1) => sec%direction=2   
-               !                   |                   Nextpoint=WestPoint(I-1,J)  => sec%direction=0   
-               !                   |==>V(I,J+1)        Nextpoint=EastPoint(I+1,J)  => sec%direction=1
-               !            U(I,J) |       U(I+1,J)
-               !            ^      |       ^
-               !   West_____|______A_______|_____EstPoint 
-               !   Point           |(I,J)        (I+1,J)  
-               !   (I-1,J)         |               
-               !                   |==>V(I,J)
-               !                   |
-               !              SoutPoint(I,J-1)
-               IF( l_oldmethod )THEN
-                  IF( zdistnew == zdistWest )      THEN
-                       sec%direction(jseg)=0 ; nextPoint = WestPoint
-                  ELSE IF( zdistnew == zdistEst )  THEN
-                       sec%direction(jseg)=1 ; nextPoint = EstPoint
-                  ELSE IF( zdistnew == zdistSouth )THEN
-                       sec%direction(jseg)=2 ; nextPoint = SouthPoint
-                  ELSE IF( zdistnew == zdistNorth )THEN
-                       sec%direction(jseg)=3 ; nextPoint= NorthPoint
-                  ENDIF
-               ELSE 
-                  IF( zdistnew == zdistWest3 )      THEN
-                       sec%direction(jseg)=0 ; nextPoint = WestPoint
-                  ELSE IF( zdistnew == zdistEst3 )  THEN
-                       sec%direction(jseg)=1 ; nextPoint = EstPoint
-                  ELSE IF( zdistnew == zdistSouth3 )THEN
-                       sec%direction(jseg)=2 ; nextPoint = SouthPoint
-                  ELSE IF( zdistnew == zdistNorth3 )THEN
-                       sec%direction(jseg)=3 ; nextPoint= NorthPoint
-                  ENDIF 
-               ENDIF
-               WRITE(cltmp,'(A11, 2(i4.4,1X) )')'nextPoint = ', nextPoint
-               CALL write_debug(jsec,cltmp)
+              !f. save next point and direction of velocity.
+              !--------------------------------------------- 
+              !nextPoint will be the one which is the closest to endingPoint.
+              !sec%direction will be direction between current point and nextPoint:
+              !It will be used to compute velocity through the segment [CurrentPoint,nextPoint}:
+              !                 
+              !A:Current Point    NorthPoint(I,J+1)   Nextpoint=NorthPoint(I,J+1) => sec%direction=3
+              !                   |                   Nextpoint=SouthPoint(I,J-1) => sec%direction=2   
+              !                   |                   Nextpoint=WestPoint(I-1,J)  => sec%direction=0   
+              !                   |==>V(I,J+1)        Nextpoint=EastPoint(I+1,J)  => sec%direction=1
+              !            U(I,J) |       U(I+1,J)
+              !            ^      |       ^
+              !   West_____|______A_______|_____EstPoint 
+              !   Point           |(I,J)        (I+1,J)  
+              !   (I-1,J)         |               
+              !                   |==>V(I,J)
+              !                   |
+              !              SoutPoint(I,J-1)
+              IF( l_oldmethod )THEN
+                 IF( zdistnew == zdistWest )      THEN
+                      sec%direction(jseg)=0 ; nextPoint = WestPoint
+                 ELSE IF( zdistnew == zdistEst )  THEN
+                      sec%direction(jseg)=1 ; nextPoint = EstPoint
+                 ELSE IF( zdistnew == zdistSouth )THEN
+                      sec%direction(jseg)=2 ; nextPoint = SouthPoint
+                 ELSE IF( zdistnew == zdistNorth )THEN
+                      sec%direction(jseg)=3 ; nextPoint= NorthPoint
+                 ENDIF
+              ELSE 
+                 IF( zdistnew == zdistWest3 )      THEN
+                      sec%direction(jseg)=0 ; nextPoint = WestPoint
+                 ELSE IF( zdistnew == zdistEst3 )  THEN
+                      sec%direction(jseg)=1 ; nextPoint = EstPoint
+                 ELSE IF( zdistnew == zdistSouth3 )THEN
+                      sec%direction(jseg)=2 ; nextPoint = SouthPoint
+                 ELSE IF( zdistnew == zdistNorth3 )THEN
+                      sec%direction(jseg)=3 ; nextPoint= NorthPoint
+                 ENDIF 
+              ENDIF
+
+              WRITE(cltmp,'(A11, 2(i4.4,1X) )')'nextPoint = ', nextPoint
+              CALL write_debug(jsec,cltmp)
         
-               !f. Save nextPoint and go to nextPoint
-               !-------------------------------------
-               prevPoint = sec%listPoint(jseg)
-               jseg = jseg+1                   !increment of number of segments that form the section
-               sec%listPoint(jseg) = nextPoint !Save next point
+              !f. Save nextPoint and go to nextPoint
+              !-------------------------------------
+              prevPoint = sec%listPoint(jseg)
+              jseg = jseg+1                   !increment of number of segments that form the section
+              sec%listPoint(jseg) = nextPoint !Save next point
    
-            ENDIF ! southP/northP/WestP/EstP == endingpoint ?
-
-         ENDDO                  !End of loop on jseg
-         sec%nb_point = jseg    !Save the number of segments that form the section
+           ENDIF ! southP/northP/WestP/EstP == endingpoint ?
+
+        ENDDO                  !End of loop on jseg
+        sec%nb_point = jseg    !Save the number of segments that form the section
 
 
      ELSE ! nb_inmesh == 0
-         DO jseg=1,nb_point_max 
-               sec%listPoint(:)=POINT_SECTION(0,0)
-         ENDDO
-         sec%direction(:)=0.
-         sec%nb_point = 0
+        DO jseg=1,nb_point_max 
+           sec%listPoint(:)=POINT_SECTION(0,0)
+        ENDDO
+        sec%direction(:)=0.
+        sec%nb_point = 0
      ENDIF
 
@@ -676,7 +660,7 @@
      CALL write_debug(jsec,"list of points in the grid : ")
      DO jseg=1,sec%nb_point 
-         ji=sec%listPoint(jseg)%I ; jj=sec%listPoint(jseg)%J
-         WRITE(cltmp, '(i4.4,X,i4.4,X,i4.4,X,f8.3,X,f8.3)' )jseg,ji,jj,glamf(ji,jj),gphif(ji,jj)
-         CALL write_debug(jsec,cltmp)
+        ji=sec%listPoint(jseg)%I ; jj=sec%listPoint(jseg)%J
+        WRITE(cltmp, '(i4.4,X,i4.4,X,i4.4,X,f8.3,X,f8.3)' )jseg,ji,jj,glamf(ji,jj),gphif(ji,jj)
+        CALL write_debug(jsec,cltmp)
      ENDDO
     
Index: /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/declarations.f90
===================================================================
--- /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/declarations.f90	(revision 2857)
+++ /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/declarations.f90	(revision 2858)
@@ -18,14 +18,9 @@
   INTEGER, PUBLIC, PARAMETER   :: numdctout=2                   ! Unit for output file
 
-  INTEGER, PUBLIC              :: jpidta,jpjdta,jpizoom,jpjzoom !dimensions des fichiers lus
-  INTEGER, PUBLIC              :: jpi,jpj,jpiglo,jpjglo         !dimensions du domaine
-  INTEGER, PUBLIC              :: nlei,nlej,nlci,nlcj
+  INTEGER, PUBLIC              :: jpi,jpj                       ! domain dimensions
   INTEGER, PUBLIC              :: nb_sec                        ! Number of section read in input file
   INTEGER, PUBLIC              :: nsecdebug = 0                 ! Number of the section to debug
 
-  INTEGER,  PUBLIC, DIMENSION(:)    , ALLOCATABLE :: mig
-  INTEGER,  PUBLIC, DIMENSION(:)    , ALLOCATABLE :: mjg  
   REAL(wp), PUBLIC ,DIMENSION(:,:)  , ALLOCATABLE :: glamf,gphif,glamt,gphit,e1t
-!  REAL(wp), PUBLIC ,DIMENSION(:,:,:), ALLOCATABLE :: tmask
   INTEGER,  PUBLIC ,DIMENSION(nb_sec_max)         :: num_sec_debug
 
@@ -54,5 +49,5 @@
                                                        zlay                ! level classes      (99 if you don't want)
      REAL(wp)                                       :: slopeSection        ! coeff directeur de la section
-     INTEGER                                        :: nb_point            ! nb de points de la section
+     INTEGER                                        :: nb_point            ! section's number of point
      TYPE(POINT_SECTION),DIMENSION(nb_point_max)    :: listPoint           ! list of point in the section
   END TYPE SECTION
Index: /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/diadct_sections.f90
===================================================================
--- /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/diadct_sections.f90	(revision 2857)
+++ /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/diadct_sections.f90	(revision 2858)
@@ -7,5 +7,5 @@
      !!
      !!
-     !! History: 2011: cbricaud Mercator-Ocean
+     !! History: 09/2011: Clement Bricaud ( Mercator-Ocean )
      !!
      !!==============================================================================
@@ -34,45 +34,4 @@
      PRINT*,'CREATION OF SECTIONS FOR NEMO diadct.F90 ROUTINE'
      PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
- 
-     !----------------------!
-     !0. Read arguments     !
-     !----------------------!
-     PRINT*,'              '
-     PRINT*,'READ ARGUMENTS'
-     PRINT*,'--------------'
- 
-     !check number of arguments and display usage message if wrong
-     narg=iargc()
-     PRINT*,'narg= ',narg
-     IF ( narg /= 4 ) THEN
-          PRINT *,' Usage : generate_sections jpidta jpjdta jpizoom jpjzoom '
-          STOP
-     ENDIF
-
-     ! read on line arguments
-     CALL getarg(1,cdum) ; READ(cdum,*) jpidta
-     CALL getarg(2,cdum) ; READ(cdum,*) jpjdta
-     CALL getarg(3,cdum) ; READ(cdum,*) jpizoom 
-     CALL getarg(4,cdum) ; READ(cdum,*) jpjzoom
-  
-     PRINT*,'jpidta  jpjdta =',jpidta,jpjdta
-     PRINT*,'jpizoom jpjzoom=',jpizoom,jpjzoom
- 
-     !------------------! 
-     !0. INITIALISATION !      
-     !------------------!     
-     PRINT*,'              '
-     PRINT*,'DOMAIN SIZE'
-     PRINT*,'--------------'
-
-     !Domain size
-     jpiglo = jpidta-jpizoom+1  ; jpjglo = jpjdta-jpjzoom+1
-     jpi    = jpiglo            ; jpj    = jpjglo
-     nlci   = jpiglo            ; nlcj   = jpjglo 
-     nlei   = jpiglo            ; nlej   = jpjglo
-
-     PRINT*,'jpiglo jpjglo = ',jpiglo,jpjglo
-     PRINT*,'jpi    jpj    = ',jpi   ,jpj 
-     PRINT*,'nlci  nlcj    = ',nlci,nlcj
  
      !-------------------!
@@ -108,4 +67,6 @@
      CALL read_coord_mesh
  
+     PRINT*,'domain sizes: '
+     PRINT*,'jpi    jpj    = ',jpi   ,jpj 
      PRINT*,'domain boundaries: '
      PRINT*,' 1   1   ',glamt(1,1),gphit(1,1)
@@ -135,7 +96,7 @@
      ENDDO
 
-     !----------------------!
-     !5.ecriture du fichier !
-     !----------------------!
+     !--------------------------------!
+     !5.Write section_ijglobal.diadct !
+     !--------------------------------!
      CALL write_sections
 
Index: /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90
===================================================================
--- /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90	(revision 2857)
+++ /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90	(revision 2858)
@@ -14,5 +14,4 @@
   
   PUBLIC  read_coord_mesh
-  PRIVATE read_ncdf
 
 CONTAINS
@@ -24,51 +23,87 @@
   !! ** Purpose :   Read a coordinate and a meshmask file in NetCDF format 
   !!---------------------------------------------------------------------      
-  !! * Local declarations
-  INTEGER                                 ::ji,jj
-
-  !!----------------------------------------------------------------------
   PRINT*,'              '
   PRINT*,'READ COORDINATES AND MESHMASK'
   PRINT*,'-----------------------------'
 
+  ! Get coordinates dimensions
+  CALL getdim(cdfile="coordinates.nc")
 
-  !Allocate coordinate and meshmask array with domain size
+  !Allocate coordinates array with domain size
   ALLOCATE(glamt(jpi,jpj)) ; ALLOCATE(gphit(jpi,jpj))
   ALLOCATE(glamf(jpi,jpj)) ; ALLOCATE(gphif(jpi,jpj))
-  ALLOCATE(e1t(jpi,jpj))   !; ALLOCATE(tmask(jpi,jpj,1))
+  ALLOCATE(e1t(jpi,jpj)  )   
 
   !Read glamt
-  CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamt",ksize=(/jpi,jpj,1,1/), &
-                   kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=glamt)
+  CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamt",ksize=(/jpi,jpj,1,1/),ptab=glamt)
 
   !Read gphit
-  CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphit",ksize=(/jpi,jpj,1,1/), &
-                   kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=gphit)
+  CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphit",ksize=(/jpi,jpj,1,1/),ptab=gphit)
 
   !Read glamf
-  CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamf",ksize=(/jpi,jpj,1,1/), &
-                   kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=glamf)
+  CALL read_ncdf(cdfile="coordinates.nc",cdvar="glamf",ksize=(/jpi,jpj,1,1/),ptab=glamf)
 
   !Read gphif
-  CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphif",ksize=(/jpi,jpj,1,1/), &
-                   kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=gphif)
+  CALL read_ncdf(cdfile="coordinates.nc",cdvar="gphif",ksize=(/jpi,jpj,1,1/),ptab=gphif)
 
   !Read e1t
-  CALL read_ncdf(cdfile="coordinates.nc",cdvar="e1t",ksize=(/jpi,jpj,1,1/), &
-                   kstart=(/jpizoom,jpjzoom,1,1/),kcount=(/jpi,jpj,1,1/),ptab=e1t)
+  CALL read_ncdf(cdfile="coordinates.nc",cdvar="e1t",ksize=(/jpi,jpj,1,1/),ptab=e1t)
 
-!  tmask(:,:,:)=1
-
-  !compute mig and mjg arrays
-  ALLOCATE(mig(jpi))
-  ALLOCATE(mjg(jpj))
-  ! local domain indices ==> data domain indices
-  DO ji = 1, jpi ; mig(ji) = ji + jpizoom - 1 ; ENDDO
-  DO jj = 1, jpj ; mjg(jj) = jj + jpjzoom - 1 ; ENDDO
-
- 
   END SUBROUTINE read_coord_mesh
 
-  SUBROUTINE read_ncdf(cdfile,cdvar,ksize,kstart,kcount,ptab)
+  SUBROUTINE getdim(cdfile)
+  !!----------------------------------------------------------------------
+  !!              ***  ROUTINE coord_mesh_read  ***
+  !!
+  !! ** Purpose :   Read a coordinate and a meshmask file in NetCDF format
+  !!----------------------------------------------------------------------
+  !! * Arguments
+  CHARACTER(*),INTENT(IN):: cdfile
+
+  !! * Local declarations
+  INTEGER           :: ncid                 ! file unit
+  INTEGER           :: idims                ! number of dimensions
+  INTEGER           :: istatus, id_var      ! dummy variable
+  CHARACTER(len=30) :: clname               ! dimension name   
+  INTEGER, ALLOCATABLE,DIMENSION(:) :: ndim ! dimensions value
+  !!----------------------------------------------------------------------
+
+  !Open file
+  istatus=NF90_OPEN(TRIM(cdfile),nf90_nowrite,ncid)
+
+  IF( istatus/= NF90_NOERR )THEN
+     PRINT*,TRIM(cdfile),' not found.stop ' ; STOP
+  ELSE
+ 
+     ! read number of dimensions   
+     istatus=NF90_INQUIRE(ncid,ndimensions=idims)
+
+     ALLOCATE( ndim(idims) )
+
+     ! read each dimension
+     PRINT*,'     File dimensions: '
+     DO id_var=1,idims
+        istatus=NF90_Inquire_Dimension(ncid,id_var,clname,ndim(id_var))
+        PRINT*,'       ',id_var,clname,ndim(id_var)
+     ENDDO
+
+     !close
+     istatus=NF90_CLOSE( ncid )              
+     IF( istatus/=NF90_NOERR )THEN 
+        PRINT*,'Problem for closing ',TRIM(cdfile);STOP
+     ELSE
+        PRINT*,'     close ',TRIM(cdfile),' OK'
+     ENDIF
+
+  ENDIF
+
+  !domain dimensions
+  jpi = ndim(1)
+  jpj = ndim(2)
+
+  DEALLOCATE( ndim )
+  END SUBROUTINE getdim
+
+  SUBROUTINE read_ncdf(cdfile,cdvar,ksize,ptab,kstart,kcount)
   !!----------------------------------------------------------------------
   !!              ***  ROUTINE coord_mesh_read  ***
@@ -81,5 +116,5 @@
   INTEGER,DIMENSION(4),INTENT(IN)                                    :: ksize
   INTEGER,DIMENSION(4),INTENT(IN),OPTIONAL                           :: kstart,kcount
-  REAL(wp),DIMENSION(ksize(1),ksize(2),ksize(3),ksize(4)),INTENT(OUT),OPTIONAL:: ptab
+  REAL(wp),DIMENSION(ksize(1),ksize(2),ksize(3),ksize(4)),INTENT(OUT):: ptab
 
   !! * Local declarations
@@ -97,11 +132,5 @@
   istatus=NF90_OPEN(TRIM(cdfile),nf90_nowrite,ncid)
   IF( istatus/= NF90_NOERR )THEN
-      PRINT*,TRIM(cdfile),' not found.stop ' ; stop
-  ELSE
-      PRINT*,'     File dimensions: '
-      DO id_var=1,4
-         istatus=NF90_Inquire_Dimension(ncid,id_var,clname,len)
-         PRINT*,'       ',id_var,clname,len
-    ENDDO
+     PRINT*,TRIM(cdfile),' not found.stop ' ; STOP
   ENDIF
 
@@ -111,5 +140,5 @@
   istatus=NF90_INQ_VARID (ncid,TRIM(cdvar),id_var)
   IF( istatus/=NF90_NOERR )THEN
-      PRINT*,TRIM(cdvar),' not found in ',TRIM(cdfile),'.stop';stop
+      PRINT*,TRIM(cdvar),' not found in ',TRIM(cdfile),'.stop';STOP
   ENDIF
 
@@ -126,5 +155,5 @@
 
   IF( istatus/=NF90_NOERR )THEN
-           PRINT*,'Problem for reading ',TRIM(cdvar),' in ',TRIM(cdfile);stop
+           PRINT*,'Problem for reading ',TRIM(cdvar),' in ',TRIM(cdfile); STOP
   ELSE
       PRINT*,'     read ',TRIM(cdvar),' OK'
Index: /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/writesections.f90
===================================================================
--- /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/writesections.f90	(revision 2857)
+++ /branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/writesections.f90	(revision 2858)
@@ -52,5 +52,4 @@
      llok=.FALSE.
      clname='section_ijglobal.diadct'
-     !CALL file_open(numdctout,clname,llok)
      CALL file_open(numdctout,clname,llok,cdform="UNFORMATTED",cdstatus="REPLACE",cdaction="WRITE")
 
@@ -77,7 +76,5 @@
              IF( secs(jsec)%nb_point .NE. 0 )THEN
                 DO jseg=1,secs(jsec)%nb_point
-                   point=POINT_SECTION( mig(secs(jsec)%listPoint(jseg)%I) , mjg(secs(jsec)%listPoint(jseg)%J) )
-                   i1 = mig(secs(jsec)%listPoint(jseg)%I) ; i2 = mjg(secs(jsec)%listPoint(jseg)%J)
-                   !WRITE(numdctout)point
+                   i1 = secs(jsec)%listPoint(jseg)%I  ; i2 = secs(jsec)%listPoint(jseg)%J
                    WRITE(numdctout)i1,i2
                 ENDDO
@@ -85,29 +82,4 @@
              ENDIF
 
-
-             !---------------------
-             !WRITE(500,*)jsec
-             !WRITE(500,*)secs(jsec)%name
-             !WRITE(500,*)secs(jsec)%llstrpond
-             !WRITE(500,*)secs(jsec)%ll_ice_section
-             !WRITE(500,*)secs(jsec)%ll_date_line
-             !WRITE(500,*)secs(jsec)%coordSec
-             !WRITE(500,*)secs(jsec)%nb_class
-             !WRITE(500,*)secs(jsec)%zsigi
-             !WRITE(500,*)secs(jsec)%zsigp
-             !WRITE(500,*)secs(jsec)%zsal
-             !WRITE(500,*)secs(jsec)%ztem
-             !WRITE(500,*)secs(jsec)%zlay
-             !WRITE(500,*)secs(jsec)%slopeSection
-             !WRITE(500,*)secs(jsec)%nb_point
-             !IF( secs(jsec)%nb_point .NE. 0 )THEN
-             !   DO jseg=1,secs(jsec)%nb_point
-             !      point=POINT_SECTION( mig(secs(jsec)%listPoint(jseg)%I) , mjg(secs(jsec)%listPoint(jseg)%J) )
-             !      WRITE(500,*)point
-             !   ENDDO
-             !   WRITE(500,*)secs(jsec)%direction(1:secs(jsec)%nb_point)
-             !ENDIF
-
-             !---------------------
           ENDDO !end of loop on sections
          
