Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/job.ksh
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/job.ksh	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/job.ksh	(revision 3104)
@@ -0,0 +1,46 @@
+#!/bin/ksh
+set -x
+#####################################################
+#     CONFIGURATION CHOICE                          #
+#####################################################
+#CONFIG=ORCA2_LIM GYRE
+CONFIG=ORCA2_LIM
+
+#####################################################
+#     CONFIGURATION DESCRIPTION                     #
+#####################################################
+# each configuration is defined by:
+#   - INSPACE: directory the coordinate netcdf file is
+#   - COOR_FIL: configuration's coordinate netcdf file
+#   - LIST:     configuration's list of sections
+#
+#####################################################
+case $CONFIG in
+#
+ORCA2_LIM)
+INSPACE=/dataref/rd/INITIALISATION/ORCA2
+COOR_FIL=coordinates.nc
+LIST=list_sections.ascii_global
+;;
+#
+GYRE)
+INSPACE=/perm/ms/fr/ar5
+COOR_FIL=mesh_mask.nc
+LIST=list_sections.ascii_GYRE
+;;
+#
+esac
+
+#####################################################
+#     RUN                                           #
+#####################################################
+BIN=diadct_sections.exe
+
+export CTLDIR=`pwd`
+cd ${CTLDIR}
+rm -f coordinates.nc ${BIN} list_sections.ascii
+cp ../${BIN} .
+ln -sf ${LIST} list_sections.ascii
+ln -s $INSPACE/$COOR_FIL ./coordinates.nc
+./${BIN}
+mv section_ijglobal.diadct section_ijglobal.diadct_$CONFIG 
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_GYRE
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_GYRE	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_GYRE	(revision 3104)
@@ -0,0 +1,3 @@
+ -80.     14.    -80.     30.   00 okstrpond noice  SEC1
+ -70.     14.    -70.     30.   00 okstrpond noice  SEC2
+ -90.     20.    -60.     20.   00 okstrpond noice  SEC3
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_exp
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_exp	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_exp	(revision 3104)
@@ -0,0 +1,27 @@
+(F7.2,1X,F7.2,1X,F7.2,1X,F7.2,1X,I2,1X,A9,1X,A5,1X,A40)
+
+long1   lat1    long2   lat2    nbclasse S/T transp? ice transp? section name
+
+-xxx.xx_-xxx.xx_-xxx.xx_-xxx.xx_ii_aaaaaaaaa_aaaaa_aaaaa....
+ 100.    -23.9   115.    -23.9  00 nostrpond noice IND_W1
+ 100.    -33.8   117.    -33.8  00 nostrpond noice IND_W2
+ 100.    -33.8   100.      2.   00 nostrpond noice IND_W3
+ 122.    -45.    122.    -30.   00 nostrpond noice IND_S1
+ 139.5   -45.    139.5   -32.   00 nostrpond noice IND_S2
+ 122.    -45.    139.5   -45.   00 nostrpond noice IND_S3
+ 150.    -28.5   165.    -28.5  00 nostrpond noice IND_E1
+ 148.    -34.    165.    -34.   00 nostrpond noice IND_E2
+ 146.    -41.7   165.    -41.7  00 nostrpond noice IND_E3
+ 165.    -41.7   165.    -28.5  00 nostrpond noice IND_E4
+ 149.    -10.    165.    -10.   00 nostrpond noice IND_E6
+ 165.    -41.7   173.    -41.7  00 nostrpond noice IND_E7
+ 115.     24.    122.     17.   00 nostrpond noice IND_N1
+ 125.      8.    134.     -2.   00 nostrpond noice IND_N2
+ 114.1    -8.45  118.9    -8.45 00 nostrpond noice IND_N3
+ 121.8    -8.7   124.1    -9.7  00 nostrpond noice IND_N4
+ 124.1    -9.7   127.    -15.   00 nostrpond noice IND_N5
+ 115.    -25.    104.     -4.   00 nostrpond noice IND_N6
+ 116.      2.3   134.     -2.   00 nostrpond noice IND_N8
+  33.   -22.2     45.    -22.2  00 nostrpond noice IND_AG1
+  25.   -33.65    44.    -33.65 00 nostrpond noice IND_AG2
+  44.   -22.2     44.    -33.65 00 nostrpond noice IND_AG3
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_global
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_global	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/list_sections.ascii_global	(revision 3104)
@@ -0,0 +1,357 @@
+ -68.    -54.5   -60.    -64.7  00 okstrpond noice ACC_Drake_Passage
+  27.    -73.     27.    -30.   00 nostrpond noice ACC_LeCap
+ 116.    -70.    116.    -33.8  00 nostrpond noice ACC_WAUS
+ 147.3   -70.    147.3   -41.   00 okstrpond noice ACC_TAS
+ 167.4   -73.    167.4   -45.5  00 nostrpond noice ACC_NZ
+ -82.     73.4   -82.     74.9  00 nostrpond   ice ARC_Lancaster_sound
+ -81.     75.1   -81.     76.6  00 nostrpond   ice ARC_Jones_Strait
+ -76.     78.4   -72.     78.3  00 nostrpond   ice ARC_Robeson_Channel
+ -64.6    60.3   -65.     61.5  00 nostrpond   ice ARC_Hudson_Strait
+  20.9    71.5    16.     77.   03 nostrpond   ice ARC_Bear_Island
+zsal
+ 0.00
+34.9
+50.0
+  33.5    68.     33.5    74.   03 nostrpond   ice ARC_Kola_Section
+zsal
+ 0.00
+34.9
+50.0
+  25.     80.     48.2    80.2  00 nostrpond   ice ARC_Spitzberg_FJLand
+  57.     80.3    65.9    76.4  03 nostrpond   ice ARC_FJLand_NovZemlija
+zsigi
+ 0.
+27.8
+30.0
+  57.     70.7    59.     70.2  00 nostrpond   ice ARC_Kara_Gate
+  25.     70.     22.6    77.8  05 nostrpond   ice ARC_Barents_Sea
+zsigi
+ 0.0
+25.0
+27.8
+29.0
+30.0
+ -21.5    79.5    15.     79.   05 nostrpond   ice ARC_Fram_Strait
+zsigi
+0.0
+25.0
+27.8
+29.0
+30.0
+ -81.     27.1   -73.     27.   04 nostrpond noice ATL_West_Atlantic_27N
+zlay
+    0.
+ 1000.
+ 6000.
+20000.
+ -35.     15.1   -35.     -8.   00 nostrpond noice ATL_Tropical_Atlantic_35W
+ -63.      9.    -55.     20.   03 nostrpond noice ATL_South_West_Atlantic
+ztem
+-2.0
+ 4.5
+40.0
+ -35.     54.    -30.5    52.   03 nostrpond noice ATL_Gibbs
+zlay
+    0.
+ 3000.
+20000.
+ -40.     -4.    -40.     10.   03 nostrpond noice ATL_North_Brazil_Current_40W
+zsigi
+ 0.0
+27.45
+30.0
+ -55.     -1.    -50.     10.   05 nostrpond noice ATL_North_Brazil_Current_10N
+zsigi
+ 0.0
+24.5
+27.125
+27.45
+30.0 
+ -53.    -30.1   -40.    -30.   00 nostrpond noice ATL_Brazil_Current_30S
+ -63     -40.1   -53.    -40.   00 nostrpond noice ATL_Falklands_Current_40S
+ -23.    -15.    -23.     15.   05 nostrpond noice ATL_Atlantic_Equatorial
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+   5.    -15.1    13.    -15.   00 nostrpond noice ATL_Benguela_Current_15S
+  10.    -40.     22.    -31.   00 nostrpond noice ATL_Agulhas_zone_30S
+  19.    -30.1   -54.0   -30.0  00 nostrpond noice ATL_Trans_Atlantic_30S
+  -7.     62.2    -4.5    58.   05 nostrpond noice ATL_England_Faroe
+zsigi
+ 0.0
+25.
+27.8
+29.
+30.0
+ -15.     65.     -7.     62.2  05 nostrpond noice ATL_Iceland_Faroe
+zsigi
+ 0.0
+25.
+27.8
+29.
+30.0
+ -37.     66.1   -22.5    66.   05 nostrpond   ice ATL_Denmark_Strait
+zsigi
+ 0.0
+25.
+27.8
+29.
+30.0
+ -44.     63.    -35.     61.   05 nostrpond   ice ATL_East_Greenland
+zsigi
+ 0.0
+25.
+27.8
+29.
+30.0
+ -45.     60.9   -45.     57.   05 nostrpond   ice ATL_South_Greenland
+zsigi
+ 0.0
+25.
+27.8
+29.
+30.0
+ -50.1    55.1   -61.     55.   05 nostrpond   ice ATL_Labrador_55N 
+zsigi
+ 0.0
+25.
+27.8
+29.
+30.0
+   8.     58.5    -2.5    57.   00 nostrpond noice ATL_England_Norway 
+  -5.7    35.5    -5.7    36.5  00 nostrpond noice ATL_Gibraltar 
+  -8.5    32.     -8.5    38.2  00 nostrpond noice ATL_Gulf_Cadiz
+ -44.     60.     -27.    56.   00 nostrpond noice ATL_Ovide1
+ -27.     56.     -16.    40.   00 nostrpond noice ATL_Ovide2
+ -16.     40.1     -8.1   40.   00 nostrpond noice ATL_Portugal
+  -4.     48.2     -8.1   43.1  05 nostrpond noice ATL_Biscay
+zsal
+ 0.0
+35.6
+36.0
+36.0
+50.0
+ -14.     26.1    -22.    26.   00 nostrpond noice ATL_South_Canaries
+ -24.     32.     -24.    37.9  03 nostrpond noice ATL_Acores_24W
+zlay
+   0.
+  800.
+20000.
+ -66.4    18.1    -63.    10.   05 nostrpond noice ATL_Antilles
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -60.1    18.1    -66.5   18.2  05 nostrpond noice ATL_Porto_Rico
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -66.5    18.1    -68.9   18.5  05 nostrpond noice ATL_Hispagnola_Porto_Rico
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -75.     20.2    -72.6   19.7  05 nostrpond noice ATL_Windward_Passage
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -77.5    18.2    -84.    14.   05 nostrpond noice ATL_Jamaica_Ridge
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -77.5    18.     -76.5   20.3  05 nostrpond noice ATL_Cuba_Jamaica
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -87.3    21.     -84.    22.2  05 nostrpond noice ATL_Yucatan_Strait_KANEC
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -80.5    22.5    -80.5   25.5  05 nostrpond noice ATL_Cuba_Florida
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -78.5    26.7    -80.5   27.   05 okstrpond noice ATL_Florida_Bahamas
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -70.     30.1    -82.    30.   05 nostrpond noice ATL_West_Atlantic_30N
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -70.     45.     -70.    30.1  05 nostrpond noice ATL_Gulf_Stream
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -55.     30.     -55.    47.9  05 nostrpond noice ATL_Gulf_Stream_2
+ztem
+-2.0
+ 4.5
+ 7.0
+12.0
+40.0
+ -55.     48.1     -3.9   48.   00 okstrpond noice ATL_Atlantic_48N
+ -85.     24.1    -10.8   24.   00 okstrpond noice ATL_Atlantic_24N
+ -75.6    24.     -10.8   24.   00 okstrpond noice ATL_Atlantic_24N1
+ -85.     26.1    -10.8   26.   04 okstrpond noice ATL_Atlantic_26N
+zlay
+    0.
+ 1000.
+ 6000.
+20000.
+ 100.    -23.9   115.    -23.9  00 nostrpond noice IND_W1
+ 100.    -33.8   117.    -33.8  00 nostrpond noice IND_W2
+ 100.    -33.8   100.      2.   00 nostrpond noice IND_W3
+ 122.    -45.    122.    -30.   00 nostrpond noice IND_S1
+ 139.5   -45.    139.5   -32.   00 nostrpond noice IND_S2
+ 122.    -45.    139.5   -45.   00 nostrpond noice IND_S3
+ 150.    -28.5   165.    -28.5  00 nostrpond noice IND_E1
+ 148.    -34.    165.    -34.   00 nostrpond noice IND_E2
+ 146.    -41.7   165.    -41.7  00 nostrpond noice IND_E3
+ 165.    -41.7   165.    -28.5  00 nostrpond noice IND_E4
+ 165.    -10.    165.    -28.5  00 nostrpond noice IND_E5
+ 149.    -10.    165.    -10.   00 nostrpond noice IND_E6
+ 165.    -41.7   173.    -41.7  00 nostrpond noice IND_E7
+ 115.     24.    122.     17.   00 nostrpond noice IND_N1
+ 125.      8.    134.     -2.   00 nostrpond noice IND_N2
+ 114.1    -8.45  118.9    -8.45 00 nostrpond noice IND_N3
+ 121.8    -8.7   124.1    -9.7  00 nostrpond noice IND_N4
+ 124.1    -9.7   127.    -15.   00 nostrpond noice IND_N5
+ 115.    -25.    104.     -4.   00 nostrpond noice IND_N6
+ 116.      2.3   134.     -2.   00 nostrpond noice IND_N8
+  33.   -22.2     45.    -22.2  00 nostrpond noice IND_AG1
+  25.   -33.65    44.    -33.65 00 nostrpond noice IND_AG2
+  44.   -22.2     44.    -33.65 00 nostrpond noice IND_AG3
+ 114.    -8.5    120.    -21.   05 nostrpond noice IND_Australia_Bali
+zlay
+    0.
+  200.
+  500.
+ 2000.
+20000.
+  45.     6.1     53.      6.   00 nostrpond noice IND_Somali_Current
+  39.   -16.      45.    -18.   00 nostrpond noice IND_Mozambique_Current
+  27.   -32.1     40.    -32.   06 nostrpond noice IND_Agulhas_Current_32S
+zsigi
+ 0.
+26.5
+26.75
+27.
+27.4
+50.0 
+  43.25   12.65   43.75   12.9  00 nostrpond noice IND_Bab_El_Manded
+  41.      0.1    50.      0.   00 nostrpond noice IND_Indian_Equatorial
+  81.    -20.     81.      8.   00 nostrpond noice IND_Indian_80E
+ 115.    -30.1   109.    -30.   00 nostrpond noice IND_West_Australia
+  52.    -24.     47.    -22.   00 nostrpond noice IND_Madagascar
+  25.    -30.1   120.    -30.   00 nostrpond noice IND_Indian_30S
+-171.5    66.2  -166.0    65.7  00 okstrpond   ice ARC_Bering_Strait
+-120.     47.1   142.5    47.   00 nostrpond noice PAC_Subarctic_Gyre_47N
+-130.     43.   -122.     47.   00 nostrpond noice PAC_OffCalif1
+-122.     40.   -130.     36.   00 nostrpond noice PAC_OffCalif2
+-118.     35.   -126.     31.   00 nostrpond noice PAC_OffCalif3
+ 121.     24.   -110.6    24.   03 okstrpond noice PAC_North_Pacific_24N
+zsigi
+20.
+27.
+30.
+-137.      7.   -137.     25.   00 nostrpond noice PAC_North_Pacific_137W_NEC
+-137.     25.   -137.     37.   00 nostrpond noice PAC_North_Pacific_137W_Kuroshio
+ 132.8    33.    137.     26.   00 nostrpond noice PAC_Kuroshio
+ 175.    -20.    175.     10.   00 nostrpond noice PAC_Pacific_Equatorial_175E
+-125.    -10.   -125.     10.   00 nostrpond noice PAC_Pacific_Equatorial_125W
+ -70.    -17.1   -76.    -17.   00 nostrpond noice PAC_Peru_Current
+ 152.    -30.1   175.    -30.   00 nostrpond noice PAC_East_Australia
+ 142.5    -8.    142.5   -13.   00 nostrpond noice PAC_Torres_Strait
+ 147.     10.    147.     50.   06 nostrpond noice PAC_Pacific_147E_10N_50N
+zsigi
+20.
+26.6
+26.7
+26.8
+27.2
+30.0
+ 165.      0.    165.     50.   06 nostrpond noice PAC_Pacific_165E
+zsigi
+20.
+26.6
+26.7
+26.8
+27.2
+30.0
+ 179.9     0.    179.9    50.   06 nostrpond noice PAC_Pacific_180E
+zsigi
+20.
+26.6
+26.7
+26.8
+27.2
+30.0
+-140.      0.   -140.     30.   06 nostrpond noice PAC_Pacific_140W
+zsigi
+20.
+26.6
+26.7
+26.8
+27.2
+30.0
+-150.     40.   -150.     60.5  00 nostrpond noice PAC_Alaskan_Gyre
+ 161.     55.   -162.     55.5  00 nostrpond noice PAC_Bering_Sea
+ 137.      0.    137.     36.   00 nostrpond noice PAC_JMA_137E_0N_35N
+ 150.    -30.    -70.    -30.   00 nostrpond noice PAC_Pacific_30S
+   0.     38.75    0.     35.36 00 nostrpond noice MED_0E
+   0.     38.75    1.37   39.   00 nostrpond noice MED_Ibiza_Channel
+   9.     36.8     9.     39.37 00 nostrpond noice MED_Sardinia_Channel
+   9.3    42.5    11.68   42.5  00 nostrpond noice MED_Corsica_Channel
+  10.5    36.6    12.9    37.85 00 nostrpond noice MED_Sicily_Strait
+  18.3    40.1    20.07   40.   00 nostrpond noice MED_Otranto_Strait
+  25.     31.5    25.     35.05 00 nostrpond noice MED_Cretan_Passage
+  28.     30.8    28.     37.32 00 nostrpond noice MED_Rhodes_Gyre
+  25.5    35.1    28.6    37.07 00 nostrpond noice MED_Kassos_Strait
+  33.4    35.     33.4    36.36 00 nostrpond noice MED_Cilician_Channel
+  22.9    36.8    23.85   35.42 00 nostrpond noice MED_Kytheron_Strait
+  10.08   60.     10.08,  56.   00 nostrpond noice BAL_WesternBaltic
+  10.     57.     13.     57.   00 nostrpond noice BAL_Kattegat
+  14.     57.     14.     53.   00 nostrpond noice BAL_Skagerrak
+ -53.5    47.     -43.    47.   00 nostrpond noice ATL_FlemishCap
+ -40.    -11.     20.    -11.1  03 okstrpond noice ATL_11S
+zsigi
+0.0
+28.12
+60.0
+  25.    -32.    120.    -32.   00 okstrpond noice IND_Indian_32S
+ 150.    -32.    -70.    -32.   00 okstrpond noice PAC_Pacific_32S
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/namelist
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/namelist	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/run/namelist	(revision 3104)
@@ -0,0 +1,8 @@
+&namdct
+!----------------------------------------!
+!nsecdebug =-1 (debug all section        !
+!          = 0 (no section to debug)     !
+!          = n (debug section # n        !
+!----------------------------------------!
+nsecdebug= 1
+/
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/compute_sections.f90	(revision 3104)
@@ -0,0 +1,692 @@
+MODULE compute_sections
+   !!=====================================================================
+   !!                       ***  MODULE  diadct  ***
+   !! Ocean diagnostics: Compute the transport through a section
+   !!
+   !! History: 2011: Clement Bricaud, Mercator-Ocean
+   !!
+   !!===============================================================
+   !! * Modules used
+   USE declarations 
+   USE sections_tools
+ 
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC compsec
+   
+   !! * Module variables
+ 
+CONTAINS
+
+  SUBROUTINE compsec(jsec,sec,lkdebug)         
+     !!---------------------------------------------------------------------
+     !!                     ***  ROUTINE compsec  ***
+     !!
+     !!  ** Purpose : Compute the series of mesh points that represent the section
+     !!               defined by its extremities.
+     !!
+     !!  ** Method  : 
+     !!          I.   Find which cells of the mesh the section is crossing
+     !!          II.  Classification of the intersections mesh/section
+     !!                  -first  classification west to east          
+     !!                  -second classification south to north 
+     !!          III. Find extremities of section in the mesh
+     !!          IV.  Find the series of mesh points that form the section
+     !!  ** Input: sec : the section to compute
+     !!
+     !!  ** Output:
+     !!---------------------------------------------------------------------
+     !! * Arguments
+     INTEGER,INTENT(IN)           :: jsec    ! number of the section
+     TYPE(SECTION), INTENT(INOUT) :: sec     ! information about the section
+     LOGICAL ,INTENT(IN)          :: lkdebug ! debug or not debug this section
+
+     !! * Local variables
+     INTEGER :: &
+        ji,jj            ,     & ! dummy loop argument
+        jseg             ,     & ! loop on segments that form the section    
+        nb_inmesh        ,     & ! number of intersections between the section and the mesh
+        nmesh                    ! number of cells in processor domain
+     INTEGER :: itest , jtest    ! dummy integer
+     REAL(wp)      :: &
+        zdistEst   , zdistNorth , zdistWest , zdistSouth ,  &! temporary scalars
+        zdistEst2  , zdistNorth2, zdistWest2, zdistSouth2,  &! temporary scalars
+        zdistEst3  , zdistNorth3, zdistWest3, zdistSouth3,  &! temporary scalars
+        zdistFirst , zdistLast  , zdistref  ,               &! temporary scalars
+        zdistante  , zdistante2 , zdistnew  , zdistnew2  ,  &! temporary scalars
+        zdeltai    , zdeltaj                                 ! temporary scalars
+     LOGICAL :: & 
+        ll_overlap_sec_left = .FALSE. , ll_overlap_sec_right = .FALSE. ,&! temporary logical
+        ll_date_domain_left = .FALSE. , ll_date_domain_right = .FALSE. ,&! temporary logical
+        ll_overlap_sec      = .FALSE. , ll_date_domain       = .FALSE. ,&! temporary logical
+        ll_test             = .FALSE.                                    ! temporary logical
+     LOGICAL :: lest, lwest, lnorth, lsouth
+     LOGICAL :: l_oldmethod
+     CHARACTER(len=120) :: cltmp
+     TYPE(COORD_SECTION) :: &
+        coord_a   , coord_b  , coord_d ,  coord_t ,               &!temporary point (long/lat)
+        coordFirst, coordLast, coordTemp                           !        "
+     TYPE(COORD_SECTION), DIMENSION(nb_point_max) :: coordSec      !intersections mesh/section    
+     TYPE(POINT_SECTION) :: &
+        endingPoint, prevPoint , nextPoint,           &!temporary point (I/J)
+        SouthPoint , NorthPoint, EstPoint , WestPoint  !        "
+     !!---------------------------------------------------------------------
+     IF( jsec==1 )THEN
+          PRINT*,'                '
+          PRINT*,'COMPUTE SECTIONS'
+          PRINT*,'----------------'
+     ENDIF
+
+     !debug
+     CALL write_debug(jsec,'compute section:')
+     CALL write_debug(jsec,'================')
+
+
+     !==================!
+     !0. Initializations!
+     !==================! 
+     
+     nmesh       = jpi*jpj          ! number of cells in processor domain
+     nb_inmesh   = 0                ! initialize number of intersections between section and the mesh
+     zdistEst  =0. ; zdistNorth=0. ; zdistWest=0. ; zdistSouth=0.   ! temporary scalars
+     zdistFirst=0. ; zdistLast =0.                                  ! temporary scalars
+     zdistante =0. ; zdistante2=0. ; zdistnew=0.  ; zdistnew2=0.    ! temporary scalars
+     zdeltai=0.    ; zdeltaj=0. 
+     coord_a     = COORD_SECTION( 0., 0. ) ; coord_b    = COORD_SECTION( 0., 0. )
+     coord_d     = COORD_SECTION( 0., 0. ) ; coord_t    = COORD_SECTION( 0., 0. ) 
+     coordFirst  = COORD_SECTION( 0., 0. ) ; coordLast  = COORD_SECTION( 0., 0. )
+     coordTemp   = COORD_SECTION( 0., 0. ) ; coordSec   = COORD_SECTION( 0., 0. )
+     endingPoint = POINT_SECTION( -1, -1 ) 
+     prevPoint   = POINT_SECTION( -1, -1 ) ; nextPoint   = POINT_SECTION( -1, -1 )
+     SouthPoint  = POINT_SECTION( -1, -1 ) ; NorthPoint  = POINT_SECTION( -1, -1 )
+     EstPoint    = POINT_SECTION( -1, -1 ) ; WestPoint   = POINT_SECTION( -1, -1 )
+
+     !===========================================================!    
+     !I. Find which cells of the mesh the section is  crossing  !
+     !===========================================================!
+
+     !loop on the mesh 
+     DO jj=2,jpj
+        DO ji=2,jpi
+           !-----------------------------------------------------------
+           !For each cell of the mesh, we find the intersection between 
+           !the section and the cell as described below:
+           !The cell is defined by 4 points A B C D 
+           !and the section by S1 and S2 
+           !
+           !Section\    ________
+           !        \ B|        |
+           !         \ |        |
+           !          \|one cell|
+           !      .....\        |
+           !           |\       |
+           !           | \      | 
+           !          A|__\_____|D
+           !       ........\
+           !----------------\-------------------------------------------
+
+           !definition of points A,B,D (we don't need C) 
+           coord_a = COORD_SECTION( glamf(ji-1,jj-1) ,gphif(ji-1,jj-1) )
+           coord_b = COORD_SECTION( glamf(ji-1,jj)   ,gphif(ji-1,jj)   )
+           coord_d = COORD_SECTION( glamf(ji,jj-1)   ,gphif(ji,jj-1)   )
+ 
+           !size of the cell 
+           zdeltai= glamf(ji-1,jj) - glamf(ji-1,jj-1) !zdeltai=[AB]
+           zdeltaj= gphif(ji,jj-1) - gphif(ji-1,jj-1) !zdeltaj=[AD]
+           
+           !intersection section/[AB]=?
+           coordTemp     = intersec(sec,coord_a,coord_b)     !compute intersection 
+ 
+           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.
+                 IF( coordTemp%lon .GE. glamf(jpi,1) .AND. coordTemp%lon .LE. glamf(jpi,1)+3) & 
+                      & ll_overlap_sec_right = .TRUE. 
+              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.
+ 
+        ENDDO
+     ENDDO !End of the loop on the mesh
+
+ 
+     !Crossing section/overlapping band (we need to know it for later):
+     !----------------------------------------------------------------- 
+     !If there is one intersection mesh/section just on the left of
+     !the overlapping band (ll_overlap_sec_left  = .TRUE.)
+     !AND there is one just the right of the overlapping band 
+     !(ll_overlap_sec_right  = .TRUE.),
+     !so the section crosses the overlapping band.
+     ll_overlap_sec = ll_overlap_sec_left .AND. ll_overlap_sec_right 
+
+     !Crossing of the domain and the date line (we need to know it for later): 
+     !------------------------------------------------------------------------
+     !If there is one point of the domain that is just on the left of the date line
+     !(ll_date_domain_left = .TRUE.) AND one point that is just on the right of the
+     !date line (ll_date_domain_right = .TRUE. ) 
+     !So the domain crosses the date line:
+     ll_date_domain = ll_date_domain_left .AND. ll_date_domain_right
+
+     !=====================================================!
+     ! II. Classification of the intersections mesh/section!
+     !=====================================================!
+
+     !   -first classification  west to east          
+     !   -second classification south to north        
+     !CAUTION: routine qcksrt doesn't work in the same way if the section
+     !and the domain crosse the date line (sec%ll_date_line=T and ll_date_domain=T)
+
+     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
+
+     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
+
+     ENDIF
+
+     !debug
+     WRITE(cltmp,'(A20,i3.3)')'number intersections = ',nb_inmesh ; CALL write_debug(jsec,cltmp)
+     CALL write_debug(jsec,'List of intersections between grid and section: ')
+     DO jseg=1,nb_inmesh
+        WRITE(cltmp,'(i4.4,1X,2(f8.3,1X) )')jseg,coordSec(jseg) ;  CALL write_debug(jsec,cltmp)
+     ENDDO
+ 
+     !=====================================================!
+     ! III. Find extremities of section in the mesh        !
+     !=====================================================!
+     !we can find section's extremities in the mesh only if
+     !there is intersection between section and mesh (nb_inmesh .ne. 0)
+
+     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,jpj
+           DO ji=1,jpi  
+              coord_t=COORD_SECTION(glamf(ji,jj),gphif(ji,jj))
+              zdistFirst = distance2(coord_t,coordFirst)
+              zdistLast = distance2(coord_t,coordLast)
+              IF( zdistFirst .LT. zdistante )THEN
+                 sec%listPoint(1) = POINT_SECTION(ji,jj)
+                 zdistante=zdistFirst
+              ENDIF
+              IF( 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
+     ENDIF 
+
+     !debug
+     CALL write_debug(jsec,"extremities of section in the grid : ")
+     ji=sec%listPoint(1)%I ; jj=sec%listPoint(1)%J
+     IF( sec%nb_point .ne. 0 )THEN
+        ji=sec%listPoint(1)%I ; jj=sec%listPoint(1)%J
+        WRITE(cltmp,'(A15,X,i4.4,X,i4.4,X,f8.3,X,f8.3)')'First point: ',sec%listPoint(1),glamf(ji,jj),gphif(ji,jj) 
+        CALL write_debug(jsec,cltmp)
+        ji=endingPoint%I ; jj=endingPoint%J
+        WRITE(cltmp,'(A15,X,i4.4,X,i4.4,X,f8.3,X,f8.3)')'Last  point: ',endingPoint,glamf(ji,jj),gphif(ji,jj)
+        CALL write_debug(jsec,cltmp)
+        !
+        coord_a=pointToCoordF(sec%listPoint(1)) ; coord_b=pointToCoordF(endingPoint)
+        ll_test = .FALSE.
+        IF ( ll_date_domain .AND. ABS( coord_a%lon - coord_b%lon ).GT. 180) ll_test= .TRUE.
+        zdistante=distance2(coord_a,coord_b ,ll_test )
+        WRITE(cltmp,'(A20,f10.3)' )'distance between IJ-extremities : ',zdistante
+        CALL write_debug(jsec,cltmp)
+        !
+        CALL write_debug(jsec,"Initial extremities : ") 
+        WRITE(cltmp,'( 2(f9.3),A3,2(f9.3) )')coordFirst,'---',coordLast
+        CALL write_debug(jsec,cltmp)
+        ll_test = .FALSE.
+        IF( ll_date_domain .AND. ABS(coordFirst%lon - coordLast%lon).GT. 180)ll_test= .TRUE.
+        zdistante=distance2(coordFirst,coordLast,ll_test)
+        WRITE(cltmp,'(A30,f10.3)')' distance between initial extremities : ',zdistante
+        CALL write_debug(jsec,cltmp)
+        CALL write_debug(jsec,"                  ")
+     ELSE
+        WRITE(cltmp,'(A50)' )"no intersection between section and mesh"
+     ENDIF
+
+     !==========================================================!
+     ! IV. Find the series of mesh points that form the section !
+     !==========================================================!
+     CALL write_debug(jsec,"Find the serie of mesh's points that form the section")
+
+     IF( sec%nb_point .ne. 0 )THEN
+
+        !The series of mesh 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 first 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  = distance2(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,102)EstPoint,WestPoint,NorthPoint,SouthPoint
+           CALL write_debug(jsec,cltmp)
+           itest=MIN(MAX(EstPoint%I,0),jpi+1) ; jtest=MIN(MAX(EstPoint%J,0),jpj+1)
+           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN 
+              WRITE(cltmp,101)'Est',glamf(itest,jtest),gphif(itest,jtest)
+              CALL write_debug(jsec,cltmp)
+           ELSE
+              CALL write_debug(jsec,"Est point out of domain")
+           ENDIF
+           !
+           itest=MIN(MAX(WestPoint%I,0),jpi+1) ; jtest=MIN(MAX(WestPoint%J,0),jpj+1)
+           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN 
+              WRITE(cltmp,101)'West',glamf(itest,jtest),gphif(itest,jtest)
+              CALL write_debug(jsec,cltmp)
+           ELSE
+              CALL write_debug(jsec,"West point out of domain")
+           ENDIF
+           !
+           itest=MIN(MAX(NorthPoint%I,0),jpi+1) ; jtest=MIN(MAX(NorthPoint%J,0),jpj+1)
+           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN 
+              WRITE(cltmp,101)'North',glamf(itest,jtest),gphif(itest,jtest)
+              CALL write_debug(jsec,cltmp)
+           ELSE
+              CALL write_debug(jsec,"North point out of domain")
+           ENDIF
+           !
+           itest=MIN(MAX(SouthPoint%I,0),jpi+1) ; jtest=MIN(MAX(SouthPoint%J,0),jpj+1)
+           IF( itest .NE. 0 .AND. itest .NE. jpi+1 .AND. jtest .NE. 0 .AND. jtest .NE. jpj+1 )THEN 
+              WRITE(cltmp,101)'South',glamf(itest,jtest),gphif(itest,jtest)
+              CALL write_debug(jsec,cltmp)
+           ELSE
+              CALL write_debug(jsec,"South point out of domain")
+           ENDIF
+           !
+           !
+100 FORMAT ( A15,2(i4.4," "),2(f8.3," ") )
+101 FORMAT ( A6,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 )
+
+               
+           !Either we are at an 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
+           !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 CAREFUL! 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 one of these lines crosses the date line, we have to add 360° 
+              ! to the longitudes of the extremities to compute the distance through
+              ! the date line and not around the Earth.
+       
+              ! 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 that 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 EastPoint 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)
+                   
+                 !be careful! 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 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
+
+
+     ELSE ! isec%nb_point == 0
+        DO jseg=1,nb_point_max 
+           sec%listPoint(:)=POINT_SECTION(0,0)
+        ENDDO
+        sec%direction(:)=0.
+        sec%nb_point = 0
+     ENDIF
+
+     !debug     
+     IF( sec%nb_point .ne. 0 )THEN
+        CALL write_debug(jsec,"-------------------------------------")
+        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)
+        ENDDO
+    
+        !test if we are an end-point
+        IF( sec%listPoint(sec%nb_point)%I .NE. endingPoint%J .AND. sec%listPoint(sec%nb_point)%J .NE. endingPoint%J )THEN   
+           PRINT*,TRIM(sec%name)," NOT ARRIVED TO endingPoint FOR jsec =  ",jsec
+        ENDIF
+     ENDIF
+
+     !now compute new slopeSection with ij-coordinates of first and last point 
+     sec%slopeSection = 0  ! default value
+     IF( sec%nb_point .ne. 0 )THEN
+        IF (  sec%listPoint(sec%nb_point)%I .NE.  sec%listPoint(1)%I ) THEN
+           sec%slopeSection = ( sec%listPoint(sec%nb_point)%J - sec%listPoint(1)%J ) /  &
+                              ( sec%listPoint(sec%nb_point)%I - sec%listPoint(1)%I )      
+        ELSE
+           sec%slopeSection = 10000._wp
+        ENDIF
+     ENDIF
+ 
+  END SUBROUTINE compsec
+
+END MODULE compute_sections 
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/declarations.f90
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/declarations.f90	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/declarations.f90	(revision 3104)
@@ -0,0 +1,61 @@
+MODULE declarations 
+  !!----------------------------------------------------------------------------
+  !!             *** declarations global varibles 
+  !!
+  !! History: 2011: Clement Bricaud, Mercator-Ocean
+  !!
+  !!----------------------------------------------------------------------------
+
+  IMPLICIT NONE
+  PUBLIC
+
+  !! * Shared module variables 
+  INTEGER, PUBLIC, PARAMETER   :: dp=8 , sp=4, wp=dp
+  INTEGER, PUBLIC, PARAMETER   :: nb_class_max  = 10            ! Max number of classes
+  INTEGER, PUBLIC, PARAMETER   :: nb_sec_max    = 150           ! Max number of sections
+  INTEGER, PUBLIC, PARAMETER   :: nb_point_max  = 2000          ! Max number of segments per section
+  INTEGER, PUBLIC, PARAMETER   :: nb_type_class = 14            ! Max number of types of classes
+  INTEGER, PUBLIC, PARAMETER   :: numnam=3                      ! Unit for namelist
+  INTEGER, PUBLIC, PARAMETER   :: numdctin=1                    ! Unit for input file
+  INTEGER, PUBLIC, PARAMETER   :: numdctout=2                   ! Unit for output file
+
+  INTEGER, PUBLIC              :: jpi,jpj                       ! domain dimensions
+  INTEGER, PUBLIC              :: nb_sec                        ! Number of sections read from input file
+  INTEGER, PUBLIC              :: nsecdebug = 0                 ! Number of the section to debug
+
+  REAL(wp), PUBLIC ,DIMENSION(:,:)  , ALLOCATABLE :: glamf,gphif,glamt,gphit,e1t
+  INTEGER,  PUBLIC ,DIMENSION(nb_sec_max)         :: num_sec_debug
+
+  TYPE POINT_SECTION
+     INTEGER :: I,J
+  END TYPE POINT_SECTION
+
+  TYPE COORD_SECTION
+     REAL(wp) :: lon,lat
+  END TYPE COORD_SECTION
+
+  TYPE SECTION
+     CHARACTER(len=60)                              :: name                ! name of the sec
+     LOGICAL                                        :: llstrpond           ! true if you want the computation of salt 
+                                                                           ! and heat transport
+     LOGICAL                                        :: ll_ice_section      ! ice surface and icevolume transport computation
+     LOGICAL                                        :: ll_date_line        ! = T if the section crosses the date-line
+     TYPE(COORD_SECTION), DIMENSION(2)              :: coordSec            ! longitude and latitude of the extremities of the section
+     INTEGER                                        :: nb_class            ! number of boundaries for density classes
+     INTEGER, DIMENSION(nb_point_max)               :: direction           ! vector direction of the point in the section
+     CHARACTER(len=40),DIMENSION(nb_class_max)      :: classname           ! characteristics of the class
+     REAL(wp), DIMENSION(nb_class_max)              :: zsigi             ,&! insitu density classes    (99 if you don't want)
+                                                       zsigp             ,&! potential density classes    (99 if you don't want)
+                                                       zsal              ,&! salinity classes   (99 if you don't want)
+                                                       ztem              ,&! temperature classes(99 if you don't want)
+                                                       zlay                ! depth level classes      (99 if you don't want)
+     REAL(wp)                                       :: slopeSection        ! slope of the section 
+     INTEGER                                        :: nb_point            ! number of points in the section
+     TYPE(POINT_SECTION),DIMENSION(nb_point_max)    :: listPoint           ! list of points in the section
+  END TYPE SECTION
+
+  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
+
+
+
+END MODULE declarations
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/diadct_sections.f90
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/diadct_sections.f90	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/diadct_sections.f90	(revision 3104)
@@ -0,0 +1,114 @@
+PROGRAM generate_sections
+     !!==============================================================================
+     !!                       ***  PROGRAM generate_sections    ***
+     !!
+     !! create a binary file containig the IJ positions of sections in global 
+     !! coordinates for the diagnostic routine diadct.F90 of NEMO
+     !!
+     !!
+     !! History: 09/2011: Clement Bricaud ( Mercator-Ocean )
+     !!
+     !!==============================================================================
+     !! * Modules used
+     USE declarations
+     USE sections_tools
+     USE readcoordmesh 
+     USE readsections
+     USE compute_sections
+     USE writesections    
+
+     IMPLICIT NONE
+ 
+     !! * Module Variables used
+     INTEGER            :: iargc, narg
+     CHARACTER(LEN=80)  :: cdum
+     INTEGER            :: jsec ,&! loop on sections
+                           jseg   ! loop on segments (parts of the section)
+     CHARACTER(len=40)  :: clname
+     LOGICAL            :: llok
+
+     NAMELIST/namdct/nsecdebug
+     !!==============================================================================     
+ 
+     PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+     PRINT*,'CREATION OF SECTIONS FOR NEMO diadct.F90 ROUTINE'
+     PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+ 
+     !-------------------!
+     !1. Read namelist   !
+     !-------------------!
+     PRINT*,'              '
+     PRINT*,'READ NAMELIST'
+     PRINT*,'--------------'
+
+     !!open, read and close namelist
+     nsecdebug=0
+     clname='namelist'
+     CALL file_open(numnam,clname,llok,cdform="FORMATTED",cdstatus="OLD",cdaction="READ")
+     IF ( llok ) THEN
+          REWIND( numnam )
+          READ  ( numnam, namdct )
+          PRINT*,'  '
+          PRINT*,'read namelist'
+          IF( nsecdebug==-1      )THEN ; PRINT*,' Debug all sections'
+          ELSE IF ( nsecdebug==0 )THEN ; PRINT*,' No section to debug'
+          ELSE IF ( nsecdebug .GE. 1 .AND. nsecdebug .LE. nb_sec_max )THEN
+              PRINT*,' Debug section number ',nsecdebug
+          ELSE
+              PRINT*,'Wrong number for nsecdebug = ',nsecdebug
+          ENDIF
+     ENDIF
+     CLOSE(numnam)
+     PRINT*,'read namelist ok'
+
+     !-------------------------------------!
+     !2. Read coordinates and meshmask     !
+     !-------------------------------------!
+     CALL read_coord_mesh
+ 
+     PRINT*,'domain sizes: '
+     PRINT*,'jpi    jpj    = ',jpi   ,jpj 
+     PRINT*,'domain boundaries: '
+     PRINT*,' 1   1   ',glamt(1,1),gphit(1,1)
+     PRINT*,' 1   jpj ',glamt(1,jpj),gphit(1,jpj)
+     PRINT*,' jpi 1   ',glamt(jpi,1),gphit(jpi,1)
+     PRINT*,'jpi jpj  ',glamt(jpi,jpj),gphit(jpi,jpj)
+
+
+
+     !----------------------!
+     !3. Read list_sections !
+     !----------------------!
+     num_sec_debug(:)=0     ! Unit numbers for debug files
+     CALL read_list_sections
+
+     !----------------------!
+     !4.Compute sections    !
+     !----------------------!
+     DO jsec=1,nb_sec
+          !we use compsec to generate the series of grid points making the section
+          IF(jsec == nsecdebug .OR. nsecdebug ==-1)THEN
+             CALL compsec(jsec,secs(jsec),.true.)
+          ELSE
+             CALL compsec(jsec,secs(jsec),.false.)
+          ENDIF
+          IF (jsec == nb_sec)PRINT*,'compute section ok '
+     ENDDO
+
+     !--------------------------------!
+     !5.Write section_ijglobal.diadct !
+     !--------------------------------!
+     CALL write_sections
+
+     !----------------------! 
+     !END                   !
+     !----------------------!
+    
+     !close debug files
+     DO jsec=1,nb_sec
+        IF( num_sec_debug(jsec) .NE. 0 )CLOSE(num_sec_debug(jsec))
+     ENDDO
+
+     PRINT*,'END END END END END END END END END END END END'
+
+END PROGRAM generate_sections
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readcoordmesh.f90	(revision 3104)
@@ -0,0 +1,199 @@
+MODULE readcoordmesh
+  !!=====================================================================
+  !!                       ***  MODULE  readcoordmesh  ***
+  !! 
+  !! History: 2011: Clement Bricaud, Mercator-Ocean
+  !! 
+  !!=====================================================================
+  !! * Modules used
+  USE netcdf
+  USE declarations
+
+  IMPLICIT NONE
+  PRIVATE
+  
+  PUBLIC  read_coord_mesh
+
+CONTAINS
+
+  SUBROUTINE read_coord_mesh
+  !!---------------------------------------------------------------------
+  !!              ***  ROUTINE coord_mesh_read  ***
+  !!
+  !! ** Purpose :   Read a coordinate file and a meshmask file in NetCDF format 
+  !!
+  !!---------------------------------------------------------------------      
+  PRINT*,'              '
+  PRINT*,'READ COORDINATES AND MESHMASK'
+  PRINT*,'-----------------------------'
+
+  ! Get coordinates dimensions
+  CALL getdim(cdfile="coordinates.nc")
+
+  !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)  )   
+
+  !Read 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/),ptab=gphit)
+
+  !Read 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/),ptab=gphif)
+
+  !Read e1t
+  CALL read_ncdf(cdfile="coordinates.nc",cdvar="e1t",ksize=(/jpi,jpj,1,1/),ptab=e1t)
+
+  END SUBROUTINE read_coord_mesh
+
+  SUBROUTINE getdim(cdfile)
+  !!----------------------------------------------------------------------
+  !!              ***  ROUTINE getdim  ***
+  !!
+  !! ** Purpose :   get dimensions of a netcdf file
+  !!
+  !!----------------------------------------------------------------------
+  !! * 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 ! dimension 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  ***
+  !!
+  !! ** Purpose :   Read a coordinate and a meshmask file in NetCDF format
+  !!
+  !!----------------------------------------------------------------------
+  !! * Arguments
+  CHARACTER(*),        INTENT(IN)                                    :: cdfile
+  CHARACTER(*),        INTENT(IN)                                    :: cdvar
+  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):: ptab
+
+  !! * Local declarations
+  INTEGER                                 ::istatus,ncid,id_var,len
+  CHARACTER(len=30) :: clname , cdvar2
+
+  INTEGER :: idims
+  INTEGER,DIMENSION(3)::idimids 
+  !!----------------------------------------------------------------------
+  ptab=0.
+  PRINT*,'read ',TRIM(cdvar),' in ',TRIM(cdfile)
+  
+  !OPEN
+  !----
+  istatus=NF90_OPEN(TRIM(cdfile),nf90_nowrite,ncid)
+  IF( istatus/= NF90_NOERR )THEN
+     PRINT*,TRIM(cdfile),' not found.stop ' ; STOP
+  ENDIF
+
+  !READ
+  !----
+  !search variable
+  istatus=NF90_INQ_VARID (ncid,TRIM(cdvar),id_var)
+  IF( istatus/=NF90_NOERR )THEN
+      PRINT*,TRIM(cdvar),' not found in ',TRIM(cdfile),'.stop';STOP
+  ENDIF
+
+  !get variable
+  !------------
+  istatus=nf90_inquire_variable(ncid,id_var, cdvar2, ndims=idims, dimids=idimids)
+  IF ( PRESENT(kstart) .AND. PRESENT(kcount) )THEN  
+      istatus=NF90_GET_VAR(ncid,id_var,ptab,start=kstart,count=kcount)
+  ELSE
+      istatus=NF90_GET_VAR(ncid,id_var,ptab)
+  ENDIF
+
+  CALL ERR_HDL(istatus)
+
+  IF( istatus/=NF90_NOERR )THEN
+           PRINT*,'Problem for reading ',TRIM(cdvar),' in ',TRIM(cdfile); STOP
+  ELSE
+      PRINT*,'     read ',TRIM(cdvar),' OK'
+  ENDIF
+
+  !CLOSE
+  !-----
+  istatus=NF90_CLOSE( ncid )
+  IF( istatus/=NF90_NOERR )THEN
+      PRINT*,'Problem for closing ',TRIM(cdfile);stop
+  ELSE
+      PRINT*,'     close ',TRIM(cdfile),' OK'
+  ENDIF
+
+
+  END SUBROUTINE read_ncdf
+
+  SUBROUTINE ERR_HDL(kstatus)
+  !! ----------------------------------------------------------
+  !!   ***  SUBROUTINE err_hdl
+  !!
+  !!   ** Purpose :  Error handle for NetCDF routine.
+  !!          Stop if kstatus indicates error conditions.
+  !!
+  !! History :
+  !!     Original: J.M. Molines (01/99)
+  !!
+  !! -----------------------------------------------------------
+  INTEGER, INTENT(in) ::  kstatus
+
+  !! -----------------------------------------------------------
+  IF( kstatus /=  NF90_NOERR ) THEN
+     PRINT *, 'ERROR in NETCDF routine, status=',kstatus
+     PRINT *,NF90_STRERROR(kstatus)
+     STOP
+  END IF
+
+  END SUBROUTINE ERR_HDL
+
+END MODULE readcoordmesh
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readsections.f90
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readsections.f90	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/readsections.f90	(revision 3104)
@@ -0,0 +1,200 @@
+MODULE readsections 
+   !!=====================================================================
+   !!                       ***  MODULE  readsections  ***
+   !! 
+   !! History: 2011: Clement Bricaud, Mercator-Ocean
+   !! 
+   !!=====================================================================
+   !! * Modules used
+   USE declarations   
+   USE sections_tools
+ 
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC read_list_sections 
+ 
+CONTAINS
+
+  SUBROUTINE read_list_sections
+     !!---------------------------------------------------------------------
+     !!         ***  ROUTINE read_list_sections ***
+     !!
+     !! ** Purpose: read ascii file 'list_sections.ascii' that contains
+     !!             section descriptions
+     !!
+     !!---------------------------------------------------------------------
+     !! * arguments
+
+     !! * Local declarations
+     INTEGER             :: jsec !loop on section number
+     INTEGER             :: iost,ji
+     INTEGER             :: iclass , jclass
+     LOGICAL             :: llok,llstrpond,llice,lldate
+     REAL(wp)            :: plon1,plat1,plon2,plat2
+     REAL(wp)            :: zslope
+     CHARACTER(len=5)    :: clclass
+     CHARACTER(len=5)    :: cdice
+     CHARACTER(len=9)    :: cdstrpond
+     CHARACTER(len=110)  :: clname,cdsecname,cltmp
+     REAL,DIMENSION(nb_type_class)     :: zclass_value
+     TYPE(COORD_SECTION)               :: coord_point1,coord_point2,coordTemp
+     TYPE(COORD_SECTION), DIMENSION(2) :: coord_sec
+     !!---------------------------------------------------------------------
+     PRINT*,'              '
+     PRINT*,'READ list_sections'
+     PRINT*,'------------------'
+     PRINT*,'               '
+
+     nb_sec=0 !initialize number of sections read in list_sections.ascii
+
+     !open and read input file
+     clname='list_sections.ascii'
+     PRINT*,'               '
+     CALL file_open(numdctin,clname,llok,cdform="FORMATTED",cdstatus="OLD",cdaction="READ")
+
+     IF ( llok ) THEN
+          PRINT*,'list_sections.ascii open '
+           PRINT*,'nb_sec_max = ',nb_sec_max
+          PRINT*,'                         '
+     
+          DO jsec=1,nb_sec_max
+
+                 !read a line corresponding to one section
+                 READ(numdctin,'(F7.2,1X,F7.2,1X,F7.2,1X,F7.2,1X,I2,1X,A9,1X,A5,1X,A40)',iostat=iost) & 
+                 &    plon1,plat1,plon2,plat2,iclass,cdstrpond,cdice,cdsecname
+                 IF (iost /= 0) EXIT  ! end of file
+
+                 ! cdsecname: change space to underscore for cdsecname
+                 cdsecname=ADJUSTL(cdsecname)
+                 ji = SCAN(TRIM(cdsecname)," ")
+                 DO WHILE(ji .NE. 0)
+                    cdsecname(ji:ji) = "_"
+                    ji = SCAN(TRIM(cdsecname)," ")
+                 ENDDO
+
+                 !computation of heat and salt transport ?
+                 llstrpond=.FALSE. ; IF( cdstrpond .EQ. 'okstrpond' ) llstrpond=.TRUE.
+
+                 !computation of ice tranpsort ?
+                 llice=.FALSE. ; IF( cdice .EQ. 'okice' ) llice=.TRUE.
+
+                 !store coordinates of the extremities
+                 coord_point1=COORD_SECTION(plon1,plat1)
+                 coord_point2=COORD_SECTION(plon2,plat2)
+                 coord_sec=(/coord_point1,coord_point2/)
+
+                 !Extremities of the section are classed
+                 lldate=.FALSE.
+                 IF(  coord_sec(2)%lon .LT. coord_sec(1)%lon  .OR.   &
+                    ((coord_sec(2)%lon .EQ. coord_sec(1)%lon) .AND.  &
+                     (coord_sec(2)%lat .LT. coord_sec(1)%lat)) ) THEN
+                      coordTemp   =coord_sec(1)
+                      coord_sec(1)=coord_sec(2)
+                      coord_sec(2)=coordTemp
+                 ENDIF
+                 IF((coord_sec(2)%lon - coord_sec(1)%lon) .GT. 180) THEN
+                      coordTemp   =coord_sec(1)
+                      coord_sec(1)=coord_sec(2)
+                      coord_sec(2)=coordTemp
+                      lldate=.TRUE.
+                 ENDIF
+
+                 !slope of the section (equidistant cylindric projection)
+                 zslope=slope_coeff(coord_sec(1),coord_sec(2),lldate)
+
+                 !!initialise global array secs
+                 secs(jsec)%llstrpond=.FALSE.  
+                 secs(jsec)%ll_date_line=.FALSE. ; secs(jsec)%nb_class=0
+                 secs(jsec)%zsigi=99.            ; secs(jsec)%zsigp=99.
+                 secs(jsec)%zsal=99.             ; secs(jsec)%ztem=99.
+                 secs(jsec)%zlay=99.
+                 secs(jsec)%nb_point=0
+
+                 !store all information in global array secs
+                 secs(jsec)%name           = cdsecname
+                 secs(jsec)%llstrpond      = llstrpond
+                 secs(jsec)%ll_ice_section = llice
+                 secs(jsec)%coordSec       = (/ coord_sec(1) , coord_sec(2) /)
+                 secs(jsec)%slopeSection   = zslope
+                 secs(jsec)%ll_date_line   = lldate
+
+                 !debug information
+                 CALL write_debug(jsec,'Informations read in ascii file:')
+                 CALL write_debug(jsec,'--------------------------------')
+                 CALL write_debug(jsec,'section name: '//secs(jsec)%name )
+                 IF( secs(jsec)%llstrpond )THEN ; CALL write_debug(jsec,'salt/heat transport computing' )
+                 ELSE                           ; CALL write_debug(jsec,'no salt/heat transport computing' )
+                 ENDIF
+                 IF( secs(jsec)%ll_ice_section )THEN ; CALL write_debug(jsec,'Ice transport computing' )
+                 ELSE                                ; CALL write_debug(jsec,'no Ice transport computing' )
+                 ENDIF
+                 WRITE(cltmp,'(A20,2f8.3)')'Extremity 1 :',secs(jsec)%coordSec(1)
+                 CALL write_debug(jsec,cltmp)
+                 WRITE(cltmp,'(A20,2f8.3)')'Extremity 2 :',secs(jsec)%coordSec(2)
+                 CALL write_debug(jsec,cltmp)
+                 WRITE(cltmp,'(A20,f8.3)')'Slope coefficient :',secs(jsec)%slopeSection
+                 CALL write_debug(jsec,cltmp)
+                 WRITE(cltmp,'(A20,i3.3)')'number of classes : ',iclass
+                 CALL write_debug(jsec,cltmp)
+                 IF( secs(jsec)%ll_date_line ) THEN ;  CALL write_debug(jsec,'section crosses date line')
+                 ELSE                                ; CALL write_debug(jsec,'section don t crosse date line')
+                 ENDIF
+                 CALL write_debug(jsec,'                        ')
+
+                 !verify number of sections and store it
+                 IF ( iclass .GT. nb_class_max) THEN
+                     PRINT*,"WARNING:  nb_class_max needs to be greater than ", iclass ; STOP
+                 ENDIF
+                 secs(jsec)%nb_class=iclass
+
+                 !read classes
+                 IF ( iclass .NE. 0 )THEN
+
+                      !classname=zsigi/zsigp/zsal/ztem/zlay
+                      READ(numdctin,'(A5)')clclass
+                      DO jclass = 1,iclass
+                         READ(numdctin,'(F9.3)',iostat=iost) zclass_value(jclass) 
+                      ENDDO 
+                      IF      ( TRIM(clclass) .EQ. 'zsigi' )THEN
+                           secs(jsec)%zsigi(1:iclass)=zclass_value(1:iclass) 
+                      ELSE IF ( TRIM(clclass) .EQ. 'zsigp' )THEN
+                           secs(jsec)%zsigp(1:iclass)=zclass_value(1:iclass) 
+                      ELSE IF ( TRIM(clclass) .EQ. 'zsal'  )THEN
+                           secs(jsec)%zsal(1:iclass)=zclass_value(1:iclass) 
+                      ELSE IF ( TRIM(clclass) .EQ. 'ztem'  )THEN
+                           secs(jsec)%ztem(1:iclass)=zclass_value(1:iclass)
+                      ELSE IF ( TRIM(clclass) .EQ. 'zlay'  )THEN
+                           secs(jsec)%zlay(1:iclass)=zclass_value(1:iclass)
+                      ELSE
+                          PRINT*,'Wrong name of class for section/clclass: ', cdsecname,TRIM(clclass)
+                      ENDIF
+
+                      IF ( jsec==nsecdebug .OR. nsecdebug==-1)THEN
+                           PRINT*,'class type = ',clclass
+                           PRINT*,'class values = ',zclass_value(1:iclass)
+                      ENDIF
+
+                 ENDIF
+
+          ENDDO !end of loop on sections
+         
+          CLOSE(numdctin) !Close file
+          IF( jsec .EQ. nb_sec_max)THEN
+              PRINT*,'   '
+              PRINT*,' nb_sec_max is less than the number of sections written in list_sections.ascii'
+              STOP
+          ELSE 
+             nb_sec=jsec-1 !number of read sections
+             PRINT*,'   '
+             PRINT*,'Number of sections read in list_sections.ascii: ',nb_sec 
+             PRINT*,'Reading of list_sections.ascii ok'
+             PRINT*,'   '
+          ENDIF
+
+     ENDIF
+
+  END SUBROUTINE  read_list_sections
+
+END MODULE readsections 
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/sections_tools.f90
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/sections_tools.f90	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/sections_tools.f90	(revision 3104)
@@ -0,0 +1,563 @@
+MODULE sections_tools
+   !!=====================================================================
+   !!                       ***  MODULE  sections_tools  ***
+   !! 
+   !! History: 2011: Clement Bricaud, Mercator-Ocean
+   !! 
+   !!=====================================================================
+   !! * Modules used
+   USE declarations
+ 
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC pointToCoordF         ! define a point with geographical coordinates
+   PUBLIC distance2             ! compute distance between 2 points
+   PUBLIC distance3             ! compute distance between a point and a line
+   PUBLIC intersec              ! intersection between 2 lines
+   PUBLIC slope_coeff           ! slope coefficient of a line
+   PUBLIC qcksrt                ! organize a list of point
+   PUBLIC write_debug           ! write debug messages
+   PUBLIC file_open             ! open a file ( ascii, bin)
+ 
+CONTAINS
+
+  TYPE(COORD_SECTION) FUNCTION pointToCoordF(p)
+     !!---------------------------------------------------------------------
+     !!         ***  FUNCTION pointToCoordF  ***
+     !!           
+     !! ** Purpose: define a point with geographical coordinates
+     !!
+     !!---------------------------------------------------------------------
+     !! * arguments
+     TYPE(POINT_SECTION), INTENT(IN) ::  p
+
+     !!---------------------------------------------------------------------
+
+     pointToCoordF = COORD_SECTION(glamf(p%I,p%J),gphif(p%I,p%J))
+
+     RETURN
+
+  END FUNCTION pointToCoordF
+
+  REAL(wp) FUNCTION distance2(coordA,coordB,ld_date_line)
+     !!---------------------------------------------------------------------
+     !!         ***  FUNCTION distance  ***
+     !!
+     !! ** Purpose:compute distance between coordA and coordB
+     !!            We add 360 to coordB%long if the line (coordA,coordB)
+     !!            crosses the date-line.
+     !!
+     !!---------------------------------------------------------------------
+     !! * arguments
+     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB
+     LOGICAL,INTENT(IN),OPTIONAL     :: ld_date_line
+
+     !! * Local declarations 
+     REAL(wp) :: zrad, zrayon
+     REAL(wp) :: zpi = 3.141592653589793
+     REAL(wp) :: zlam1, zlam2, zphi1, zphi2
+
+     !-----------------------------------------------
+     zrayon = 6367000.0
+     zrad   = zpi/180.
+     zlam1  = coordA%lon*zrad
+     zlam2  = coordB%lon*zrad
+     zphi1  = coordA%lat*zrad
+     zphi2  = coordB%lat*zrad
+     !
+     !function output
+     distance2=0.001*2.*zrayon *ASIN(SQRT(SIN((zphi2-zphi1)/2.0)**2.0+COS(zphi1)*COS(zphi2)*SIN((zlam2-zlam1)/2.0)**2.0) )
+     !
+     RETURN
+  END FUNCTION distance2
+
+  REAL(wp) FUNCTION distance3(coordA,coordB,coordC,ll_debug) 
+     !!---------------------------------------------------------------------
+     !!         ***  FUNCTION cosinus  ***
+     !!
+     !! ** Purpose: compute distance between a point B and a line (AC)
+     !!
+     !! ** Methode: use Al-Kashi theorem
+     !!
+     !!              B                   A first point of the section
+     !!             / \                  B intermediate point on the section
+     !!            / | \                 C last point of the section
+     !!           /  |  \               
+     !!          /   |   \               angle=ACB
+     !!         /    |    \              
+     !!        /     |     \             distance3=Bd
+     !!       /      |      \
+     !!    C  -------d------- A  
+     !!
+     !!---------------------------------------------------------------------
+     !! * arguments
+     TYPE(COORD_SECTION), INTENT(IN) :: coordA,coordB,coordC
+     LOGICAL                         :: ll_debug
+
+     !! * Local declarations
+     REAL(wp)                        :: ztmp
+     REAL(wp)                        :: za, zb, zc, zangle
+
+     !-----------------------------------------------
+     za = SQRT( (coordB%lon-coordC%lon)**2 + (coordB%lat-coordC%lat)**2 )
+     zb = SQRT( (coordA%lon-coordC%lon)**2 + (coordA%lat-coordC%lat)**2 )
+     zc = SQRT( (coordA%lon-coordB%lon)**2 + (coordA%lat-coordB%lat)**2 )
+     !
+     IF( za /= 0. .AND. zb /= 0. )THEN
+        ztmp = ( za**2 + zb**2  - zc**2 ) / ( 2.0*za*zb )
+        ztmp = MIN(ztmp,1.00_wp)
+        IF( ztmp==1.00_wp )THEN ; zangle = 0.0 
+        ELSE                    ; zangle = ABS(acos(ztmp)) 
+        ENDIF
+     ELSE
+        PRINT*,'You should not have this situation: za zb =',za,zb ; STOP
+     ENDIF
+     !
+     distance3 = za*ABS(sin(zangle)) 
+     !
+     !function output
+     RETURN
+  END FUNCTION distance3
+
+  REAL(wp) FUNCTION slope_coeff(coordA,coordB,ld_dateline)
+     !!---------------------------------------------------------------------
+     !!         ***  Function slope_coeff  ***
+     !!
+     !! ** Purpose: Compute slope coefficient of the line (coordA,coordB)
+     !!
+     !! ** Method:
+     !!           Usual method :
+     !!                   slope_coeff = (latB-latA)/(lonB-lonA)
+     !!           Special case: the segment [A,B] crosses the date-line (lddate=T)
+     !!                   slope_coeff = (latB-latA)/(360+lonB-lonA)
+     !!---------------------------------------------------------------------
+     !! * arguments
+     TYPE(COORD_SECTION), INTENT(IN) :: coordA, coordB
+     LOGICAL,INTENT(IN),OPTIONAL     :: ld_dateline
+
+     !! * Local declarations
+     REAL(wp)                        :: zcoeff
+     INTEGER                         :: idateline
+
+     !!---------------------------------------------------------
+     !initialization
+     zcoeff = 1.e20
+     idateline=0
+     IF( PRESENT( ld_dateline))THEN
+       IF( ld_dateline )idateline=1
+     ENDIF
+
+     !compute slope coefficient
+     IF ( coordB%lon-coordA%lon .NE. 0.) &
+         zcoeff=(coordB%lat-coordA%lat)/(360*idateline+coordB%lon-coordA%lon)
+
+     !output
+     slope_coeff = zcoeff
+     RETURN
+
+  END FUNCTION slope_coeff 
+
+  TYPE(COORD_SECTION) FUNCTION intersec(sec,coord_a,coord_b)
+     !!---------------------------------------------------------------------
+     !!         ***  Function intersec  ***
+     !!
+     !! ** Purpose:Return the coordinates of the intersection point, 
+     !!            between segment [a,b] and the section.
+     !!            If no intersection the point = -9999.
+     !!
+     !! ** Method:(coord_a,coord_b) => y=za2*x+zb2  (2)
+     !!            sec              => y=za1*x+zb1  (1)
+     !!    Intersection = (X,Y) solves (1) and (2)
+     !!    Vefify that (X,Y) is in [coord_a,coord_b] AND in the section
+     !!         
+     !! ** Action: 1. compute za1, za1 
+     !!            2. compute zb1, zb2
+     !!            3. compute X and Y
+     !!            4.Verify that (zX,zY) is in [coord_a,coord_b] 
+     !!              and in [sec%coordSec(1),sec%coordSec(2)]    
+     !!
+     !! History: Author: 10/05 Matthieu Laborie
+     !!          Additions:
+     !!           05-2007: (C Bricaud) add special cases
+     !!                                (crossing date line)
+     !----------------------------------------------------------------------------- 
+     !! * arguments
+     TYPE(SECTION), INTENT(IN)        :: sec
+     TYPE(COORD_SECTION)  , INTENT(IN):: coord_a, coord_b
+
+     !! * Local declarations
+     TYPE(COORD_SECTION) :: coordInter
+     REAL(wp)            :: za1,za2, &! slope coefficient
+                            zb1,zb2, &! 
+                            zX,zY     ! coordinates of intersection
+     LOGICAL             :: ll_date_line=.FALSE. !segment [a,b] crosses the date-line ?
+
+     !----------------------------------------------------------------------------
+
+     !=================!
+     !0. INITALIZATION !
+     !=================!
+     coordInter=COORD_SECTION(-9999.,-9999.) 
+     ll_date_line=.FALSE.    
+     !we need to know if [coord_a,coord_b] crosses the date-line 
+     if(coord_a%lon-coord_b%lon .GT. 180)ll_date_line=.TRUE.
+
+     !=======================!
+     !1. compute za1 and za2 !
+     !=======================!
+     za1=sec%slopeSection
+     za2=slope_coeff(coord_a,coord_b,ll_date_line)
+ 
+     !=======================!
+     !2. compute zb1 and zb2 !
+     !=======================!
+
+     ! Compute coefficient b for a straight line y=a*x+b
+     !        Usual method:  knowing value of a, we compute b with coordinates of 1 point:
+     !                       b=latA-a.lonA or b=latB-a.lonB
+     !        Particular case: the straight line crosses the date line; so it is in 2 parts:
+     !                         one on the left of the date-line and one the right
+     !                         then,b can could have 2 values:
+     ! As the date-line can be crossed by the section and the segment of the mesh [a,b],
+     ! we have to check if the segment [a,b] crosses the date-line (we know it already for the
+     ! section(sec%ll_date_line): ll_date_line
+     ! Then, there are 4 cases for computing b1 and b2 ( sec%ll_date_line=T/F AND ll_date_line=T/F )
+     !=========================================================================================!
+     ! CASE A :                               !CASE C:
+     ! sec%ll_date_line=T AND ll_date_line =T !sec%ll_date_line=F AND ll_date_line =T
+     !                                        !     |  -180 !+180   |
+     ! this case doesn't exist                !     |       !       |
+     !========================================! a___|_______!_______|__b
+     !CASE B:                                 !     |       !       |
+     !sec%ll_date_line=T AND ll_date_line=F   !     |       !       |
+     !     \     |                            ! section  date-line  section
+     !   a__\__b |                            ! b1: usual method
+     !       \   |                            ! b2: depend of longitude of
+     !        \  |                            !     sec%coordSec(1) et de  sec%coordSec(2)
+     !         \ |                            !=================================================!
+     !          \|                            !CASE D:
+     !     +180  \  -180                      !sec%ll_date_line=F AND ll_date_line =F
+     !           |\                           !      |       !       |
+     !           | \                          !      | +180  ! -180  |
+     !           |  \                         !  a___|___b   !   a___|___b
+     !           |a__\__b                     !      |       !       |
+     !           |    \                       !      |       !       |
+     !    date-line   section                 !   section date-line section
+     !b2: usual method                        !b1: usual method
+     !b1: depend of longitude of a and b.     !b2: usual method
+     !==========================================================================================!
+
+     IF( sec%ll_date_line )THEN
+          IF( ll_date_line )THEN
+              !CASE A: this case doesn't exist
+              !-------
+              zb1=1.e20
+              zb2=1.e20
+          ELSE 
+              !CASE B:
+              !-------
+              !compute zb2: 
+              !zb2 = coord_b%lat - za2 * coord_b%lon = coord_a%lat - za2 * coord_a%lon 
+              zb2=coord_b%lat - za2 * coord_b%lon  
+
+              !compute zb1:
+              IF( coord_a%lon .GT. 0 .AND. coord_b%lon .GT. 0  )THEN
+                  zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
+              ELSE  
+                  zb1 = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
+              ENDIF 
+          ENDIF
+     ELSE
+          IF( ll_date_line )THEN
+              !CASE C:
+              !-------
+              !Compute zb1: 
+              !zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon 
+              !    = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon 
+              zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
+
+              !Compute zb2:
+              IF( sec%coordSec(1)%lon .GT. 0 .AND. sec%coordSec(2)%lon .GT. 0 )THEN
+                  zb2 = coord_a%lat - za2 * coord_a%lon
+              ELSE  
+                  zb2 = coord_b%lat - za2 * coord_b%lon
+              ENDIF
+          ELSE
+              !CASE D:
+              !-------
+              !compute zb1:
+              !zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon  
+              !    = sec%coordSec(2)%lat - za1 * sec%coordSec(2)%lon
+              zb1 = sec%coordSec(1)%lat - za1 * sec%coordSec(1)%lon
+
+              !compute zb2:
+              !zb2 = coord_b%lat - za2 * coord_b%lon = coord_a%lat - za2 * coord_a%lon
+              zb2 = coord_b%lat - za2 * coord_b%lon  
+          ENDIF
+     ENDIF 
+
+     IF( (za1 - za2) .NE. 0) THEN
+          !================================!
+          !3. Compute intersection (zX,zY) !
+          !================================!
+          IF( za1 == 1.e20 ) THEN ! Case X=constant
+              zX = sec%coordSec(1)%lon
+              zY = za2 * zX + zb2
+          ELSE IF ( za2 == 1.e20) THEN ! Case X=constant
+              zX = coord_b%lon
+              zY = za1 * zX + zb1
+          ELSE ! Case zY=A*zX+B
+              zX = (zb2 - zb1) / (za1 - za2)
+              zY = (za1 * zX ) + zb1
+          ENDIF
+           
+          !==============================================!
+          !4.Verify that (zX,zY) is in [coord_a,coord_b] !
+          !  and in [sec%coordSec(1),sec%coordSec(2)]    !
+          !==============================================!
+          !Be careful! The test is not the same for all configurations
+          IF( sec%ll_date_line )THEN
+
+              IF( ll_date_line )THEN
+                  !CASE A: this case doesn't exist
+                  !-------
+                  coordInter=COORD_SECTION(-9999.,-9999.)
+              ELSE
+                  !CASE B:
+                  !-------   
+                  IF( zX .GE. MIN(coord_a%lon,coord_b%lon ) .AND. &
+                      zX .LE. MAX(coord_a%lon,coord_b%lon ) .AND. &
+                      zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
+                      zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
+                    ((zX .LE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
+                      zX .GE. -180) .OR.  &
+                     (zX .GE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
+                      zX .LE. 180)) .AND. &
+                      zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
+                      zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN
+                     coordInter = COORD_SECTION(zX, zY)
+                  ENDIF
+              ENDIF
+          ELSE
+
+              IF( ll_date_line )THEN
+                  !CASE C:
+                  !-------
+                  IF( ((zX .LE. MIN(coord_a%lon,coord_b%lon ) .AND. &
+                        zX .GE. -180 )  .OR. &
+                       (zX .GE. MAX(coord_a%lon,coord_b%lon ) .AND. &
+                        zX .LE. 180 )) .AND. & 
+                        zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
+                        zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
+                        zX .GE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
+                        zX .LE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
+                        zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
+                        zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN 
+                        coordInter = COORD_SECTION(zX, zY)   
+                  ENDIF 
+              ELSE
+                  !CASE D:
+                  !-------
+                  IF( zX .GE. MIN(coord_a%lon,coord_b%lon ) .AND. &
+                      zX .LE. MAX(coord_a%lon,coord_b%lon ) .AND. &
+                      zY .GE. MIN(coord_a%lat,coord_b%lat ) .AND. &
+                      zY .LE. MAX(coord_a%lat,coord_b%lat ) .AND. &
+                      zX .GE. MIN(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
+                      zX .LE. MAX(sec%coordSec(1)%lon,sec%coordSec(2)%lon ) .AND. &
+                      zY .GE. MIN(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) .AND. &
+                      zY .LE. MAX(sec%coordSec(1)%lat,sec%coordSec(2)%lat ) ) THEN
+                      coordInter = COORD_SECTION(zX, zY)
+                  ENDIF             
+              ENDIF
+          ENDIF   
+            
+     ENDIF
+      
+     !output
+     intersec = coordInter
+     RETURN
+   
+  END FUNCTION intersec
+
+  SUBROUTINE qcksrt(arr1,arr2,n)
+     !!---------------------------------------------------------------------
+     !!         ***  SUBROUTINE qcksrt  ***
+     !!
+     !! ** Purpose : organize a list point by latitude , and after by longitude
+     !!              using "numerical recipies " routine.
+     !!---------------------------------------------------------------------
+     !! * arguments
+     INTEGER,INTENT(IN)                 :: n           ! number of points
+     REAL(wp),DIMENSION(n),INTENT(INOUT):: arr1,arr2   ! longitude and latitude
+
+     !! * Local declarations
+     INTEGER , PARAMETER                :: m = 7, nstack = 500
+     REAL(wp), PARAMETER                :: fm=7875.,fa=211.,fc=1663.,fmi=1./fm
+     INTEGER                            :: il, ir, i, iq, ist    !local integers
+     INTEGER                            :: jj                    !loop indices
+     REAL(wp)                           :: zfx, za, zb           !local real
+     INTEGER,DIMENSION(nstack)          :: istack(nstack)        !temp array
+     !---------------------------------------------------
+     ist=0
+     il=1
+     ir=n
+     zfx=0.
+     !
+ 10  IF( ir-il .LT. m )THEN
+
+        DO jj = il+1,ir
+           za = arr1(jj)
+           zb = arr2(jj)
+           DO i = jj-1,1,-1
+              IF( arr1(i) .LE. za )GOTO 12
+              arr1(i+1) = arr1(i)
+              arr2(i+1) = arr2(i)
+           ENDDO
+           i=0
+ 12        arr1(i+1) = za
+           arr2(i+1) = zb
+        ENDDO
+        IF( ist .EQ. 0 )RETURN
+        ir  = istack(ist)
+        il  = istack(ist-1)
+        ist = ist-2
+     ELSE
+        i  = il
+        jj = ir
+        zfx = MOD(zfx*fa+fc,fm)
+        iq = il+(ir-il+1)*(zfx*fmi)
+        za  = arr1(iq)
+        arr1(iq) = arr1(il)
+        zb = arr2(iq)
+        arr2(iq) = arr2(il)
+ 20     CONTINUE
+ 21     IF( jj .GT. 0 )THEN
+           IF( za .LT. arr1(jj) )THEN
+              jj = jj-1
+              GOTO 21
+           ENDIF
+        ENDIF
+        IF( jj .LE. i )THEN
+           arr1(i) = za
+           arr2(i) = zb
+           GOTO 30
+        ENDIF
+        arr1(i) = arr1(jj)
+        arr2(i) = arr2(jj)
+        i=i+1
+ 22     IF( i .LE. n )THEN
+           IF( za .GT. arr1(i) )THEN
+              i = i+1
+              GOTO 22
+           ENDIF
+        ENDIF 
+        IF( jj .LE.i )THEN
+           arr1(jj) = za
+           arr2(jj) = zb
+           i        = jj
+           GOTO 30
+        ENDIF
+        arr1(jj) = arr1(i)
+        arr2(jj) = arr2(i)
+        jj = jj - 1
+        GOTO 20
+ 30     ist = ist + 2
+        IF( ist .GT. nstack )PAUSE 'nstack too small'
+        IF( ir-i .GE. i-1 )THEN
+           istack(ist  ) = ir
+           istack(ist-1) = i + 1
+           ir = i-1
+        ELSE
+           istack(ist  ) = i-1
+           istack(ist-1) = il
+           il = i+1
+        ENDIF
+     ENDIF
+     GOTO 10
+  END SUBROUTINE qcksrt
+
+  SUBROUTINE write_debug(ksec,cd_write)
+     !!---------------------------------------------------------------------
+     !!         ***  SUBROUTINE write_debug  ***
+     !!
+     !! ** Purpose: write debug messages
+     !!
+     !!---------------------------------------------------------------------
+     !! * arguments
+     INTEGER           :: ksec       !number of section
+     CHARACTER(len=*)  :: cd_write   !message to write
+
+     !! * Local declarations
+     INTEGER           :: iunit
+     LOGICAL           :: llok
+     CHARACTER(len=80) :: clfilename
+   
+     !!---------------------------------------------------------------------
+     IF( ksec == nsecdebug .OR. nsecdebug == -1 )THEN
+
+
+        !open / verify is debug output file is open
+
+        clfilename=TRIM(secs(ksec)%name)//"_debug"
+        iunit = num_sec_debug(ksec)
+
+        IF( iunit .EQ. 0 )THEN
+           PRINT*,"Open debug file: "//TRIM(clfilename)
+           iunit = 100 + ksec
+           num_sec_debug(ksec) = iunit 
+ 
+           CALL file_open(iunit,clfilename,llok,cdform="FORMATTED",cdstatus="REPLACE",cdaction="WRITE")
+        
+           IF( .NOT. llok )THEN 
+               PRINT*,"Can not open TRIM(clfilename)." ; STOP
+           ENDIF   
+        ENDIF
+      
+       WRITE(iunit,*)TRIM(cd_write)
+
+     ENDIF
+
+  END SUBROUTINE write_debug
+
+  SUBROUTINE file_open(knum,cdfile,ldok,cdform,cdstatus,cdaction)
+     !!---------------------------------------------------------------------
+     !!         ***  ROUTINE file_open ***
+     !!
+     !! ** Purpose: open a file
+     !!
+     !!---------------------------------------------------------------------
+     !! * arguments
+     INTEGER          ,INTENT(IN)  :: knum       ! file unit number
+     CHARACTER(len=*),INTENT(IN)   :: cdfile     ! file name to open
+     LOGICAL          ,INTENT(OUT) :: ldok       ! =.TRUE. if file exists and is corectly opened
+     CHARACTER(len=*),INTENT(IN)   :: cdform,   &! FORM arguments for OPEN function
+                                      cdstatus, &! STATUS arguments for OPEN function
+                                      cdaction   ! ACTION arguments for OPEN function
+                                      
+     !! * Local declarations
+     INTEGER :: iost
+     LOGICAL :: llbon=.FALSE. !check existence of file
+     !!---------------------------------------------------------------------
+     ldok = .FALSE.                        ! initialization; file not found and not opened
+
+     !check presence of file
+     IF( cdstatus .EQ. "OLD" ) THEN
+        INQUIRE( FILE=cdfile, EXIST=llbon ) ! check presence of namelist
+        IF( llbon )THEN ;  PRINT*,TRIM(cdfile)//' EXISTS'
+        ELSE            ;  PRINT*,TRIM(cdfile)//' NOT EXISTS' ; STOP
+        ENDIF
+     ENDIF
+
+     !open file
+     OPEN(UNIT=knum,FILE=TRIM(cdfile),FORM=cdform,status=cdstatus,action=cdaction,iostat=iost)
+     IF ( iost == 0 )THEN
+          ldok=.TRUE.
+     ELSE
+          PRINT*,TRIM(cdfile)//' bad opening. STOP.'; STOP
+     ENDIF
+
+  END SUBROUTINE file_open
+
+END MODULE sections_tools 
Index: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/writesections.f90
===================================================================
--- branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/writesections.f90	(revision 3104)
+++ branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/TOOLS/SECTIONS_DIADCT/src/writesections.f90	(revision 3104)
@@ -0,0 +1,91 @@
+MODULE writesections 
+   !!=====================================================================
+   !!                       ***  MODULE  writesections  ***
+   !! 
+   !! History: 2011: Clement Bricaud, Mercator-Ocean
+   !! 
+   !!=====================================================================
+   !! * Modules used
+   USE declarations   
+   USE sections_tools
+ 
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC  write_sections 
+   PRIVATE file_open 
+ 
+CONTAINS
+
+  SUBROUTINE write_sections
+     !!---------------------------------------------------------------------
+     !!         ***  ROUTINE read_list_sections ***
+     !!
+     !! ** Purpose
+     !!
+     !! ** Method
+     !!
+     !! ** Input
+     !!
+     !! ** Action
+     !!
+     !! History
+     !!---------------------------------------------------------------------
+     !! * arguments
+
+     !! * Local declarations
+     INTEGER             :: jsec ,&!loop on sections 
+                            jseg   !loop on segments
+     INTEGER             :: i1, i2 !temporary integers
+     LOGICAL             :: llok   !
+     CHARACTER(len=40)   :: clname !
+     TYPE(POINT_SECTION) :: point  !coordinates of a point
+
+     !!---------------------------------------------------------------------
+
+     PRINT*,'              '
+     PRINT*,'WRITE SECTIONS'
+     PRINT*,'--------------'
+
+     !open output file
+     llok=.FALSE.
+     clname='section_ijglobal.diadct'
+     CALL file_open(numdctout,clname,llok,cdform="UNFORMATTED",cdstatus="REPLACE",cdaction="WRITE")
+
+     !print informations
+     IF ( llok ) THEN
+          PRINT*,TRIM(clname),' open. '
+
+          DO jsec=1,nb_sec
+  
+             WRITE(numdctout)jsec
+             WRITE(numdctout)secs(jsec)%name
+             WRITE(numdctout)secs(jsec)%llstrpond
+             WRITE(numdctout)secs(jsec)%ll_ice_section
+             WRITE(numdctout)secs(jsec)%ll_date_line
+             WRITE(numdctout)secs(jsec)%coordSec
+             WRITE(numdctout)secs(jsec)%nb_class
+             WRITE(numdctout)secs(jsec)%zsigi
+             WRITE(numdctout)secs(jsec)%zsigp
+             WRITE(numdctout)secs(jsec)%zsal
+             WRITE(numdctout)secs(jsec)%ztem
+             WRITE(numdctout)secs(jsec)%zlay
+             WRITE(numdctout)secs(jsec)%slopeSection
+             WRITE(numdctout)secs(jsec)%nb_point
+             IF( secs(jsec)%nb_point .NE. 0 )THEN
+                DO jseg=1,secs(jsec)%nb_point
+                   i1 = secs(jsec)%listPoint(jseg)%I  ; i2 = secs(jsec)%listPoint(jseg)%J
+                   WRITE(numdctout)i1,i2
+                ENDDO
+                WRITE(numdctout)secs(jsec)%direction(1:secs(jsec)%nb_point)
+             ENDIF
+
+          ENDDO !end of loop on sections
+         
+          CLOSE(numdctout) !Close file
+     ENDIF
+
+  END SUBROUTINE  write_sections
+
+END MODULE writesections 
