diff --git a/tools/check_commit_ial.sh b/tools/check_commit_ial.sh
index 2bc421e6dad9cb8c210819102c0f6deb8a33a273..771c7015fd40b8fb219e42be0962788c99c02b4c 100755
--- a/tools/check_commit_ial.sh
+++ b/tools/check_commit_ial.sh
@@ -721,10 +721,14 @@ if [ $check -eq 1 ]; then
         fi
       fi
       if [ $t -eq 0 ]; then
+        cmd="$PHYEXTOOLSDIR/compare.py"
+        if [ ! -x $cmd ]; then
+          echo "Command not found: \"$cmd\""
+          exit 10
+        fi
         if [ $(basename $file) == ICMSHFPOS+0002:00 ]; then
           #historic files
-          cmd="cmp $file1 $file2 256 256"
-          output='stderr'
+          cmd="$cmd --binary $file1 $file2 256"
         elif [ $(basename $file) == DHFDLFPOS+0002 ]; then
           #DDH files
           ddh_images="$HOMEPACK/$name/ddh_diff_${tag}.png"
@@ -732,32 +736,15 @@ if [ $check -eq 1 ]; then
             [ ! -d /d0/images/$USER ] && mkdir /d0/images/$USER
             ddh_images="$ddh_images /d0/images/$USER/ddh_diff_${tag}.png"
           fi
-          cmd="$PHYEXTOOLSDIR/comp_DDH.py"
-          if [ ! -x $cmd ]; then
-            echo "Command not found: \"$cmd\""
-            exit 10
-          fi
-          cmd="$cmd $file1 $file2 $ddh_images"
-          output='stdout'
+          cmd="$cmd --ddh $file1 $file2 --ddhplots $ddh_images"
         elif [ $(basename $file) == NODE.001_01 ]; then
           #Output listing
-          cmd="$PHYEXTOOLSDIR/diffNODE.001_01"
-          if [ ! -x $cmd ]; then
-            echo "Command not found: \"$cmd\""
-            exit 11
-          fi
-          cmd="$cmd $file1 $file2 --norm-max-diff=0."
-          output='stdout'
+          cmd="$cmd --node $file1 $file2"
         else
-          cmd="cmp $file1 $file2"
-          output='stderr'
+          cmd="$cmd --binary $file1 $file2 0"
         fi
         set +e
-        if [ $output == 'stderr' ]; then
-            mess=$($cmd 2>&1)
-        else
-            mess=$($cmd 2>/dev/null)
-        fi
+        mess=$($cmd 2>/dev/null)
         t=$?
         set -e
       fi
diff --git a/tools/check_commit_mesonh.sh b/tools/check_commit_mesonh.sh
index e310d9694332d294e5451e9466a746668de0c8d4..b690ef9d10b7edf2f497b93091ea585445453f31 100755
--- a/tools/check_commit_mesonh.sh
+++ b/tools/check_commit_mesonh.sh
@@ -518,26 +518,25 @@ if [ $check -eq 1 ]; then
       #Comparison
       if [ -f $file1 -a -f $file2 ]; then
         # Compare variable of both Synchronous and Diachronic files with printing difference
-        echo "Compare with python..."
+        echo "Comparison for case $t..."
         set +e
         if [ "$file3" == "" ]; then
-          $PHYEXTOOLSDIR/compare.py --f1 $file1 --f2 $file2
-          t=$?
+          $PHYEXTOOLSDIR/compare.py --backup $file1 $file2
+          r=$?
         else
-          $PHYEXTOOLSDIR/compare.py --f1 $file1 --f2 $file2 --f3 $file3 --f4 $file4
-          t=$?
+          $PHYEXTOOLSDIR/compare.py --backup $file1 $file2 --diac $file3 $file4
+          r=$?
         fi
         set -e
-        allt=$(($allt+$t))
+        allt=$(($allt+$r))
 
         #Check bit-repro before date of creation of Synchronous file from ncdump of all values
         #(pb with direct .nc file checks)
-        echo "Compare with ncdump..."
         set +e
-        diff <(ncdump $file1 | head -c $bit_diff) <(ncdump $file2 | head -c $bit_diff)
-        t=$?
+        $PHYEXTOOLSDIR/compare.py --ncdump $file1 $file2 $bit_diff
+        r=$?
         set -e
-        allt=$(($allt+$t))
+        allt=$(($allt+$r))
       else
         [ ! -f $file1 ] && echo "  $file1 is missing"
         [ ! -f $file2 ] && echo "  $file2 is missing"
diff --git a/tools/check_commit_testprogs.sh b/tools/check_commit_testprogs.sh
index 6a0ba9aa489ca3d43873585228fa88f822c8011f..938e1ca65c063eff14bd5a3e0dd0ce711db13062 100755
--- a/tools/check_commit_testprogs.sh
+++ b/tools/check_commit_testprogs.sh
@@ -385,13 +385,9 @@ if [ $check -eq 1 ]; then
       fi
       if [ $te -eq 0 ]; then
         set +e
-        mess=$(cmp <(cat $file1 | sed 's/\.\.//g' | sed 's/~=//g' | sed 's/!=//g' | grep -v 'Total time: ' | sed 's/-0.00000E+00|/ 0.00000E+00|/g' | sed 's/-0.00000E+00 / 0.00000E+00 /g' | sed 's/-0.00000E+00-/ 0.00000E+00-/g') \
-                   <(cat $file2 | sed 's/\.\.//g' | sed 's/~=//g' | sed 's/!=//g' | grep -v 'Total time: ' | sed 's/-0.00000E+00|/ 0.00000E+00|/g' | sed 's/-0.00000E+00 / 0.00000E+00 /g' | sed 's/-0.00000E+00-/ 0.00000E+00-/g') 246 246 2>&1)
+        mess=$($PHYEXTOOLSDIR/compare.py --testprogs $file1 $file2)
         te=$?
         set -e
-        #The use of "<()" bash syntax replaces the actual file name seen by cmp
-        #We modify the cmp output to display the actual file names
-        mess=$(echo $mess | sed "s#^.*differ# $file1 $file2 differ#")
       fi
       [ $te -ne 0 ] && message="$message $mess \n"
       alltests=$(($alltests+$te))
diff --git a/tools/comp_DDH.py b/tools/comp_DDH.py
deleted file mode 100755
index a2a625fe18e5a038357f3edaebe06c3bff1b0085..0000000000000000000000000000000000000000
--- a/tools/comp_DDH.py
+++ /dev/null
@@ -1,93 +0,0 @@
-#!/usr/bin/env python3
-
-import matplotlib
-matplotlib.use('Agg')
-import os
-import shutil
-import epygram
-import numpy
-import matplotlib.pyplot as plt
-epygram.init_env()
-
-def comp_DDH(filename1, filename2, output_fig, tol_ad=3E-7, tol_rd=1.E-6, verbose=False):
-    if filename1 != filename2:
-        r1 = epygram.formats.resource(filename1, 'r')
-        r2 = epygram.formats.resource(filename2, 'r')
-    
-        l1 = set(r1.listfields())
-        l2 = set(r2.listfields())
-        pb_var = len(l1.symmetric_difference(l2)) != 0
-    
-        def comp(fid, v1, v2):
-            t = numpy.all(v1 == v2)
-            toplt = False
-            if not t:
-                if verbose: print(fid, ':')
-                if numpy.array(v1).ndim == 0:
-                    v1 = numpy.array([v1])
-                    v2 = numpy.array([v2])
-                for i in range(len(v1)):
-                    if v1[i] - v2[i] != 0.:
-                        ad = v1[i] - v2[i]
-                        rd = 200 * (v1[i] - v2[i]) / (v1[i] + v2[i])
-                        if verbose: print("   v1={v1}, v2={v2}, diff={ad}, rdiff={rd}".format(v1=v1[i], v2=v2[i], ad=ad, rd=rd))
-                        if abs(ad) > tol_ad and abs(rd) > tol_rd:
-                          if verbose: print("  ==> plot")
-                          toplt = True
-            return fid if toplt else None
-        toplt = []
-        for fid in [fid for fid in l1.intersection(l2) if fid != 'DOCFICHIER']:
-            v1 = r1.readfield(fid)
-            v2 = r2.readfield(fid)
-            if isinstance(v1, epygram.base.FieldSet):
-                for i in range(len(v1)): #fieldset
-                    toplt.append(comp(fid, v1[i].getdata(), v2[i].getdata()))
-            else:
-                toplt.append(comp(fid, v1.getdata(), v2.getdata()))
-        toplt = [fid for fid in toplt if fid is not None]
-        pb_val = len(toplt) > 0
-        if pb_val:
-            ncols, nrows = min(10, len(toplt)), 1+(len(toplt)-1)//10
-            figure, ax = plt.subplots(ncols=ncols, nrows=nrows,
-                                      figsize=(5 * ncols, 10 * nrows), squeeze=False)
-            ax = ax.flatten()
-            figure.suptitle(filename1 + ' ' + filename2)
-            for ifid, fid in enumerate(toplt):
-                v1 = r1.readfield(fid)
-                v2 = r2.readfield(fid)
-                assert(len(v1) == len(v2))
-                for i in range(len(v1)): #fieldset
-                    ad = v1[i].getdata() - v2[i].getdata()
-                    ax[ifid].plot(v1[i].getdata(), v1[i].geometry.vcoordinate.levels, label='v1')
-                    ax[ifid].plot(v2[i].getdata(), v2[i].geometry.vcoordinate.levels, label='v2')
-                    ax[ifid].legend()
-                    ax[ifid].twiny().plot(ad, v1[i].geometry.vcoordinate.levels, label='diff', color='black', ls=':')
-                    ad = numpy.abs(ad)
-                    rd = (200 * numpy.abs(v1[i].getdata() - v2[i].getdata()) / numpy.abs(v1[i].getdata() + v2[i].getdata()))
-                    rd = rd[ad != 0.].max()
-                    ad = ad.max()
-                    ax[ifid].set_title("{fid}:\nmax_ad={ad}\nmax_rd={rd}%".format(fid=fid, ad=ad, rd=rd))
-            figure.savefig(output_fig[0])
-            for filename in output_fig[1:]:
-                 shutil.copyfile(output_fig[0], filename)
-        if pb_var and pb_val:    
-            message = "Variables are different and values of common variables are also different"
-        elif pb_var:
-            message = "Variables are different but values of common variables are equal"
-        elif pb_val:
-            message = "Values are different"
-        else:
-            message = ""
-        if pb_val:
-            message += ", plot is available in the folowing file(s): " + ', '.join(output_fig)
-    else: #same file
-        pb_var = False
-        pb_val = False
-        message = ""
-
-    print(message)
-    return 1 if pb_var or pb_val else 0
-
-if __name__ == '__main__':
-    import sys
-    sys.exit(comp_DDH(sys.argv[1], sys.argv[2], sys.argv[3:]))
diff --git a/tools/compare.py b/tools/compare.py
index 16410e00d3faddfd4647bb5f6a8daa7eb0e9a6f1..f210d1d14937d3bb698cb79a3f790ec43d8e903e 100755
--- a/tools/compare.py
+++ b/tools/compare.py
@@ -1,25 +1,36 @@
 #!/usr/bin/env python3
 
+import matplotlib
+matplotlib.use('Agg')
+import os
+os.environ['NUMEXPR_MAX_THREADS'] = '1'
+import shutil
+import epygram
+import numpy
+import matplotlib.pyplot as plt
+epygram.init_env()
 import xarray as xr
 import sys
+import subprocess
+import difflib
+import re
 
+#List of budgtes groups to compare
 avail_groups=['Stations/sta1',
-'LES_budgets/Miscellaneous/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/Mean/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/Resolved/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/Subgrid/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/Surface/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/BU_KE/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/BU_THL2/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/BU_WTHL/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/BU_RT2/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/BU_WRT/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'LES_budgets/BU_THLR/Cartesian/Not_time_averaged/Not_normalized/cart/',
-'Budgets/TH','Budgets/UU','Budgets/WW',
-'Budgets/RV','Budgets/RI','Budgets/RC',
-'Budgets/RG','Budgets/RS','Budgets/RH','Budgets/TK']
-
-tol_ad=1E-12 # Error max for budgets
+              'LES_budgets/Miscellaneous/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/Mean/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/Resolved/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/Subgrid/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/Surface/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/BU_KE/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/BU_THL2/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/BU_WTHL/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/BU_RT2/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/BU_WRT/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'LES_budgets/BU_THLR/Cartesian/Not_time_averaged/Not_normalized/cart/',
+              'Budgets/TH','Budgets/UU','Budgets/WW',
+              'Budgets/RV','Budgets/RI','Budgets/RC',
+              'Budgets/RG','Budgets/RS','Budgets/RH','Budgets/TK']
 
 def compareBACKUPFiles(file_user, file_ref):
   status = 0
@@ -61,7 +72,7 @@ def compareBACKUPFiles(file_user, file_ref):
       pass
   return status
 
-def compareTSERIESFiles(file_user, file_ref):
+def compareTSERIESFiles(file_user, file_ref, tol_ad=1E-12):
   status = 0
   da = xr.open_dataset(file_user)
   da2 = xr.open_dataset(file_ref)
@@ -110,22 +121,321 @@ def compareTSERIESFiles(file_user, file_ref):
       pass
   return status
 
+def comp_DDH(filename1, filename2, output_fig, tol_ad=3E-7, tol_rd=1.E-6, verbose=False):
+    if filename1 != filename2:
+        r1 = epygram.formats.resource(filename1, 'r')
+        r2 = epygram.formats.resource(filename2, 'r')
+    
+        l1 = set(r1.listfields())
+        l2 = set(r2.listfields())
+        pb_var = len(l1.symmetric_difference(l2)) != 0
+    
+        def comp(fid, v1, v2):
+            t = numpy.all(v1 == v2)
+            toplt = False
+            if not t:
+                if verbose: print(fid, ':')
+                if numpy.array(v1).ndim == 0:
+                    v1 = numpy.array([v1])
+                    v2 = numpy.array([v2])
+                for i in range(len(v1)):
+                    if v1[i] - v2[i] != 0.:
+                        ad = v1[i] - v2[i]
+                        rd = 200 * (v1[i] - v2[i]) / (v1[i] + v2[i])
+                        if verbose: print("   v1={v1}, v2={v2}, diff={ad}, rdiff={rd}".format(v1=v1[i], v2=v2[i], ad=ad, rd=rd))
+                        if abs(ad) > tol_ad and abs(rd) > tol_rd:
+                          if verbose: print("  ==> plot")
+                          toplt = True
+            return fid if toplt else None
+        toplt = []
+        for fid in [fid for fid in l1.intersection(l2) if fid != 'DOCFICHIER']:
+            v1 = r1.readfield(fid)
+            v2 = r2.readfield(fid)
+            if isinstance(v1, epygram.base.FieldSet):
+                for i in range(len(v1)): #fieldset
+                    toplt.append(comp(fid, v1[i].getdata(), v2[i].getdata()))
+            else:
+                toplt.append(comp(fid, v1.getdata(), v2.getdata()))
+        toplt = [fid for fid in toplt if fid is not None]
+        pb_val = len(toplt) > 0
+        if pb_val and output_fig is not None and len(output_fig) > 0:
+            ncols, nrows = min(10, len(toplt)), 1+(len(toplt)-1)//10
+            figure, ax = plt.subplots(ncols=ncols, nrows=nrows,
+                                      figsize=(5 * ncols, 10 * nrows), squeeze=False)
+            ax = ax.flatten()
+            figure.suptitle(filename1 + ' ' + filename2)
+            for ifid, fid in enumerate(toplt):
+                v1 = r1.readfield(fid)
+                v2 = r2.readfield(fid)
+                assert(len(v1) == len(v2))
+                for i in range(len(v1)): #fieldset
+                    ad = v1[i].getdata() - v2[i].getdata()
+                    ax[ifid].plot(v1[i].getdata(), v1[i].geometry.vcoordinate.levels, label='v1')
+                    ax[ifid].plot(v2[i].getdata(), v2[i].geometry.vcoordinate.levels, label='v2')
+                    ax[ifid].legend()
+                    ax[ifid].twiny().plot(ad, v1[i].geometry.vcoordinate.levels, label='diff', color='black', ls=':')
+                    ad = numpy.abs(ad)
+                    rd = (200 * numpy.abs(v1[i].getdata() - v2[i].getdata()) / numpy.abs(v1[i].getdata() + v2[i].getdata()))
+                    rd = rd[ad != 0.].max()
+                    ad = ad.max()
+                    ax[ifid].set_title("{fid}:\nmax_ad={ad}\nmax_rd={rd}%".format(fid=fid, ad=ad, rd=rd))
+            figure.savefig(output_fig[0])
+            for filename in output_fig[1:]:
+                 shutil.copyfile(output_fig[0], filename)
+        if pb_var and pb_val:    
+            message = "Variables are different and values of common variables are also different"
+        elif pb_var:
+            message = "Variables are different but values of common variables are equal"
+        elif pb_val:
+            message = "Values are different"
+        else:
+            message = ""
+        if pb_val and output_fig is not None and len(output_fig) > 0:
+            message += ", plot is available in the folowing file(s): " + ', '.join(output_fig)
+    else: #same file
+        pb_var = False
+        pb_val = False
+        message = ""
+
+    print(message)
+    return 1 if pb_var or pb_val else 0
+
+default_spnorms='VORTICITY,DIVERGENCE,TEMPERATURE,KINETIC ENERGY'
+default_gpnorms='HUMI.SPECIFI,CLOUD_WATER,ICE_CRYSTAL,CLOUD_FRACTI'
+def comp_NODE(f1, f2, spnorms=None, gpnorms=None, norm_max_diff=0.050000):
+    """
+    Adapted from the Ryad's diffNODE.001_01
+    :param f1, f2: filenames to compare
+    :param spnorms: coma separated list of spectral variables. '' to not
+                    compare spectral norms. None to use the default
+                    (""" + default_spnorms + """)
+    :param gpnorms: coma separated list of gridpoint variables. '' to not
+                    compare gridpoint norms. None to use the default
+                    (""" + default_gpnorms + """)
+    :param norm_max_diff: maximum difference allowed
+    """
+    
+    from collections import defaultdict
+
+    if spnorms is None: spnorms = default_spnorms
+    if gpnorms is None: gpnorms = default_gpnorms    
+    opts = {'spnorms': spnorms.split(','),
+            'gpnorms': gpnorms.split(','),
+            'norm_max_diff': norm_max_diff}
+
+    def xave(f):
+        with open(f, 'r') as fh:
+            gpregs = []
+            if opts['gpnorms']:
+                if len(opts['gpnorms']) == 1 and opts['gpnorms'][0] == '*':
+                    gpregs = [re.compile(r'^\s*GPNORM\s+\b(\S+)\b')]
+                else:
+                    gpregs = [re.compile(rf'^\s*GPNORM\s+\b({re.escape(gp)})\b') for gp in opts['gpnorms']]
+            x = []
+            for line in fh:
+                line = line.strip()
+                if re.match(r'^\s*GPNORM\s+', line):
+                    for gpreg in gpregs:
+                        match = gpreg.match(line)
+                        if match:
+                            F = match.group(1)
+                            line = next(fh)
+                            if not re.match(r'\s*AVE\s+', line):
+                                continue
+                            line = line.strip().lstrip('AVE')
+                            values = line.split()
+                            x.extend([[F, value] for value in values])
+                            break
+                elif re.match(r'^\s*SPECTRAL\s+NORMS\s+-\s+', line):
+                    spnormk = []
+                    spnormv = []
+                    while True:
+                        line = next(fh)
+                        if not re.match(r'\s+LEV\s+', line):
+                            break
+                        index = {spnorm: line.find(spnorm) for spnorm in opts['spnorms']}
+                        index = {k:v for (k, v) in index.items() if v > 0}
+                        spnormk = sorted([sp for sp in opts['spnorms'] if sp in index.keys()],
+                                         key=lambda spnorm: index[spnorm] if index[spnorm] >= 0 else float('inf'))
+                        line = next(fh)
+                        if not re.match(r'\s+AVE\s+', line):
+                            break
+                        spnormv = line.split()[1:] #'1:' to remove 'AVE'
+                        for spnormk, spnormv in zip(spnormk, spnormv):
+                            x.append([spnormk, spnormv])
+            return x
+
+    def center(s, n):
+        i = 0
+        while len(s) < n:
+            s = f" {s}" if i % 2 else f"{s} "
+            i += 1
+        return s
+
+    def main(f1, f2):
+    
+        fx1 = xave(f1)
+        fx2 = xave(f2)
+
+        title = "************* NORMS DIFFERENCES *************"
+        print(center(title, 118))
+        print(center('=' * len(title), 118))
+        print()
+    
+        x = [[]]
+        diff = {}
+        zero = 0
+        numb = 0
+    
+        tag1 = "NORMDIFF"
+        tag2 = "NORMSTAT"
+    
+        nout = 0
+        while fx1 and fx2:
+            (field1, x1) = fx1.pop(0)
+            (field2, x2) = fx2.pop(0)
+    
+            if field1 != field2:
+                return 1
+                #sys.exit("Field mismatch {} != {}".format(field1, field2))
+    
+            x1, x2 = x1.strip(), x2.strip()
+            if x1 and x2:
+                dx = float(x1) - float(x2)
+                dr = (2 * dx) / (float(x1) + float(x2)) if float(x1) + float(x2) > 0 else 0.0
+    
+                sdx = '{:17.9e}'.format(dx)
+                sdr = '{:17.9e}'.format(dr)
+    
+                dx = float(sdx)
+                dr = float(sdr)
+    
+                x[-1].append(" {} | {:20} | {:17.9e}  |  {:17.9e}  |  {:17}  |  {:17} {}\n".format(
+                    tag1, center(field1, 20), float(x1), float(x2), sdx, sdr, '*' if dr > opts['norm_max_diff'] else ''))
+    
+                nout += 1 if abs(dr) > opts['norm_max_diff'] else 0
+    
+                if abs(dr) > 0:
+                    n = int(numpy.floor((numpy.log(abs(dr)) / numpy.log(10))))
+                    diff[n] = diff.get(n, 0) + 1
+                else:
+                    zero += 1
+    
+                numb += 1
+            else:
+                x.append([])
+    
+        print(" {} |                      |{:19} | {:19} | {:19} | {:19}".format(
+            tag1, center("NORM(REF)", 19), center("NORM(EXP)", 19),
+            center("NORM(REF)-NORM(EXP)", 19), center("(NORM(REF)-NORM(EXP))", 19)))
+    
+        print(" {} |                      |{:19} | {:19} | {:19} | {:19}".format(
+            tag1, '', '', '', center("/NORM(REF)", 19)))
+    
+        for i in range(len(x)):
+            if not x[i]:
+                break
+            print("".join(x[i]))
+    
+        print("\n")
+    
+        diff_cumul = 0
+        perc_cumul = 0
+        for n1 in sorted(diff.keys()):
+            n2 = n1 + 1
+            diff_val = diff[n1]
+            perc = 100 * diff_val / numb
+            diff_cumul += diff_val
+            perc_cumul += perc
+            print(" {} | {:3d} .. {:3} | {:3d} / {:3d} | {:3d} / {:3d} | {:6.2f} %, {:6.2f} %\n".format(
+                  tag2, n1, n2, diff_val, numb, diff_cumul, numb, perc, perc_cumul))
+    
+        if nout:
+            print("\n")
+            text = "WARNING : SOME NORMS DIFFERENCES ARE OUTSIDE ALLOWED LIMIT OF {:6.2f} %\n".format(
+                100 * opts['norm_max_diff'])
+            print(text * 5)
+        return nout
+    
+    return main(f1, f2)
+
+def comp_binary(f1, f2, offset):
+    #pyhton filcmp does not allow to specify an offset
+    offset = [str(offset), str(offset)] if offset != 0 else []
+    p = subprocess.run(['cmp', f1, f2] + offset, capture_output=True, encoding='UTF8') 
+    if p.returncode != 0:
+        print(p.stdout)
+    return p.returncode
+
+def comp_ncdump(f1, f2, nbytes):
+    ncdumps = [subprocess.run(['ncdump', f], capture_output=True, encoding='UTF8').stdout[:int(nbytes)]
+               for f in (f1, f2)]
+    print(''.join(difflib.unified_diff(*[ncdumps[i].splitlines(keepends=True) for i in (0, 1)])))
+    return 0 if ncdumps[0] == ncdumps[1] else 1
+
+def comp_testprogs(f1, f2):
+    def read(f):
+        with open(f, 'r') as fd:
+            s = fd.read()
+            for p in ('\.\.', '~=', '!='): s = re.sub(p, '', s)
+            s = re.sub(r'\-0.00000E\+00([|\- ])', r' 0.00000E+00\1', s)  
+            s = re.sub(r'\n\sTotal time:.*\n', '\n', s)
+        return s[247:]
+    r = read(f1) == read(f2)
+    if not r:
+        print("{f1} {f2} differ".format(f1=f1, f2=f2))
+    return 0 if r else 1
+
 if __name__ == "__main__":
    import argparse
    import sys
-   parser = argparse.ArgumentParser(description='Compare toutes les variables si trouvées dans les fichiers backup et time series')
-   value = argparse.ArgumentParser()
-   parser.add_argument('--f1', metavar='file1', type=str, help="Backup file1 user ")
-   parser.add_argument('--f2', metavar='file2', type=str, help="Backup file2 reference")
-   parser.add_argument('--f3', metavar='file3', type=str, help=".000 file1 user ")
-   parser.add_argument('--f4', metavar='file4', type=str, help=".000 file2 reference")
+   parser = argparse.ArgumentParser(description='Compare all the variables in the different files')
+   parser.add_argument('--backup', metavar=('BACKUP_USER', 'BACKUP_REF'), nargs=2,
+                       type=str, help="Backup files (user and reference)")
+   parser.add_argument('--diac', metavar=('DIAC_USER', 'DIAC_REF'), nargs=2,
+                       type=str, help="Diachronic .000 files (user and reference)")
+   parser.add_argument('--ddh', nargs=2, metavar=('DDH_USER', 'DDH_REF'),
+                       type=str, help="DDH files (user and reference)")
+   parser.add_argument('--ddhplots', metavar='DDHPLOT', nargs='+',
+                       help="Plot filenames for DDH differences")
+   parser.add_argument('--node', nargs=2, metavar=('DDH_USER', 'DDH_REF'),
+                       type=str, help="NODE files (user and reference) to compare norms")
+   parser.add_argument('--binary', nargs=3, metavar=('BIN_USER', 'BIN_REF', 'OFFSET'),
+                       type=str, help="Binary files (user and reference) and offset")
+   parser.add_argument('--ncdump', nargs=3, metavar=('NC_USER', 'NC_REF', 'BYTES'),
+                       type=str, help="Netcdf files (user and reference) whose ncdump output" + \
+                                      "first bytes must be compared, with number of bytes")
+   parser.add_argument('--testprogs', nargs=2, metavar=('LST_USER', 'LST_REF'),
+                       type=str, help="Testprogs listing (user and ref)")
    args = parser.parse_args()
-   totalstatus=0
-   status1=compareBACKUPFiles(args.f1, args.f2)
-   totalstatus += status1
-   print('status1 = ' + str(status1))
-   if args.f3:
-     status2=compareTSERIESFiles(args.f3, args.f4)
-     totalstatus += status2
-     print('status2 = ' + str(status2))
+
+   totalstatus = 0
+   if args.backup:
+       status = compareBACKUPFiles(*args.backup)
+       totalstatus += status
+       print('status for backup files = ' + str(status))
+   if args.diac:
+       status = compareTSERIESFiles(*args.diac)
+       totalstatus += status
+       print('status for diachronic files = ' + str(status))
+   if args.ddh:
+       status = comp_DDH(*args.ddh, args.ddhplots)
+       totalstatus += status
+       #print('status for ddh files = ' + str(status))
+   if args.node:
+       status = comp_NODE(*args.node, norm_max_diff=0.)
+       totalstatus += status
+       #print('status for NODE files = ' + str(status))
+   if args.binary:
+       status = comp_binary(*args.binary)
+       totalstatus += status
+       #print('status for binary files = ' + str(status))
+   if args.ncdump:
+       status = comp_ncdump(*args.ncdump)
+       totalstatus += status
+       print('status for ncdump of files = ' + str(status))
+   if args.testprogs:
+       status = comp_testprogs(*args.testprogs)
+       totalstatus += status
+       #print('status for testprogs listings = ' + str(status))
    sys.exit(totalstatus)
diff --git a/tools/diffNODE.001_01 b/tools/diffNODE.001_01
deleted file mode 100755
index d33640104415782a37f87e7c0c49acd84be02a80..0000000000000000000000000000000000000000
--- a/tools/diffNODE.001_01
+++ /dev/null
@@ -1,428 +0,0 @@
-#!/usr/bin/perl -w
-
-use strict;
-use FindBin qw ($Bin);
-use lib $Bin;
-use FileHandle;
-#Frame was inlined in this script
-#use Tools::Frame;
-use Data::Dumper;
-use POSIX qw (floor);
-use Getopt::Long;
-
-my %opts = (
-            'spnorms' => 'VORTICITY,DIVERGENCE,TEMPERATURE,KINETIC ENERGY',
-            'gpnorms' => '',
-            'norm-max-diff' => .050000,
-            'jo-max-diff'   => .030000,
-           );
-
-sub frame
-{
-  my ($t, $width, $border) = @_;
-  my $len = length ($t);
-
-  $border = 1 unless (defined ($border));
-
-  my $S = ' ';
-
-  my ($C, $V, $H) = ('*', '|', '-');
-
-  ($C, $V, $H) = (' ', ' ', ' ')
-    unless ($border);
-
-  my $df = 3;
-  $width ||= $len + 2 * $df;
-
-  my $line1 = $C . ($H x ($width-2)) . $C;
-  my $line2 = $V . ($S x ($width-2)) . $V;
-
-
-  my $TEXT = '';
-
-  $TEXT .= "$line1\n";
-  for (1 .. ($df-1)/2)
-    {
-      $TEXT .= "$line2\n";
-    }
-
-  die ("Cannot frame text: `$t'\n")
-    if ($width - 2 * $df <= 0);
-  
-  while ($t)
-    {
-      my $s = substr ($t, 0, $width - 2 * $df, '');
-
-      my $i = 0;
-      while (length ($s) < $width - 2 * $df)
-        {
-          if ($i % 2)
-            {
-              $s = " $s";
-            }
-          else
-            {
-              $s = "$s ";
-            }
-          $i++;
-        }
-      my $linet = $V . ($S x ($df-1)) . $s .  ($S x ($df-1)) . $V;
-      $TEXT .= "$linet\n";
-    }
-
-  for (1 .. ($df-1)/2)
-    {
-      $TEXT .= "$line2\n";
-    }
-  $TEXT .= "$line1\n";
-}
-
-
-sub xave
-{
-  my $f = shift;
-  my $fh = 'FileHandle'->new ("<$f");
-
-  my @gpregs;
-
-  if ($opts{gpnorms})
-    {
-      if ((scalar (@{$opts{gpnorms}}) == 1) && ($opts{gpnorms}[0] eq '*'))
-        {
-          @gpregs = (qr/^\s*GPNORM\s+\b(\S+)\b/);
-        }
-      else
-        {
-          @gpregs = map { qr/^\s*GPNORM\s+\b(\Q$_\E)\b/ } @{$opts{gpnorms}};
-        }
-    }
-
-  my @x;
-  MAIN: while (defined (my $line = <$fh>))
-    {
-      AGAIN:
-
-      if ($line =~ m/^\s*GPNORM\s+/o)
-        {
-          for my $gpreg (@gpregs)
-            {
-              if ($line =~ $gpreg)
-                {
-                  my $F = $1;
-                  $line = <$fh>;
-                  next MAIN  unless ($line =~ m/AVE\s+0/o);
-                  for ($line)
-                    {
-                      s/^\s*AVE\s*//o; 
-                      s/\s+/\n/go; 
-                    }
-                  push @x, map { [ $F, $_ ] } split (m/\n/o, $line);
-                  next MAIN;
-                }
-            }
-        }
-
-      if ($line =~ s/^\s*SPECTRAL\s+NORMS\s+-\s+//o)
-        {
-          AGAIN_SPNORMS:
-
-
-          goto AGAIN
-            unless (($line = <$fh>) =~ s/^\s+LEV\s+//o);
-
-          my %index;
-          %index = ();
-          for my $spnorm (@{$opts{spnorms}})
-            {
-              my $index = index ($line, $spnorm);
-              $index{$spnorm} = $index 
-                if ($index >= 0);
-            }
-
-          my @spnormk = sort { $index{$a} <=> $index{$b} } 
-                        grep { defined $index{$_} } 
-                        @{$opts{spnorms}};
-
-          goto AGAIN
-            unless (($line = <$fh>) =~ s/^\s+AVE\s+//o);
-
-          my @spnormv = split (m/\s+/o, $line);
-
-          while (@spnormk)
-            {
-              my $spnormk = shift (@spnormk);
-              my $spnormv = shift (@spnormv);
-              die ("$spnormk, $spnormv\n")
-                unless (defined ($spnormk) && defined ($spnormv));
-              push @x, [ $spnormk, $spnormv ];
-            }
-
-          goto AGAIN_SPNORMS;
-
-        }
-    }
-
-  return @x;
-}
-
-sub xobstype
-{
-  my $f = shift;
-  my $fh = 'FileHandle'->new ("<$f");
-
-  # keep final value for obs number & JO
-
-  my @x;
-  while (defined (my $line = <$fh>))
-    {
-      next unless ($line =~ m/ObsType\s+(\d+)\s+Total:\s*(\d+)\s+(\S+)\s+/o);
-      @{$x[$1]}{qw (number JO)} = ($2, $3);
-    }
-
-  return @x;
-}
-
-sub xjog
-{
-  my $f = shift;
-  my $fh = 'FileHandle'->new ("<$f");
-
-  # keep final value for obs number & JO
-
-  my %x;
-  while (defined (my $line = <$fh>))
-    {
-      next unless ($line =~ m/Jo Global\s*:\s*(\d+)\s+(\S+)/o);
-      @x{qw (number JO)} = ($1, $2);
-    }
-
-  return \%x;
-}
-
-sub center
-{
-  my ($s, $n) = @_;
-  my $i = 0;
-  while (length ($s) < $n) 
-    {
-      $s = $i % 2 ? " $s" : "$s ";
-      $i++;
-    }
-  return $s;
-}
-
-
-
-&GetOptions 
-  ('spnorms=s'       => \$opts{'spnorms'},
-   'gpnorms=s'       => \$opts{'gpnorms'},
-   'norm-max-diff=s' => \$opts{'norm-max-diff'},
-   'jo-max-diff=s'   => \$opts{'jo-max-diff'},
-  );
-
-$opts{'spnorms'} = [ split (m/,/o, $opts{'spnorms'}) ];
-$opts{'gpnorms'} = [ split (m/,/o, $opts{'gpnorms'}) ];
-
-my ($f1, $f2) = @ARGV;
-
-die ("Usage: $0 NODE.001_01 NODE.001_01.ref\n")
-  unless ($f1 && $f2);
-
-my @fx1 = &xave ($f1);
-my @fx2 = &xave ($f2);
-
-print &frame ("NORMS DIFFERENCES", 121);
-
-my @x = ([]);
-my %diff;
-my $zero = 0;
-my $numb = 0;
-
-my $tag1 = "NORMDIFF";
-my $tag2 = "NORMSTAT";
-
-my $nout = 0;
-
-while (defined (my $fx1 = shift (@fx1)) && defined (my $fx2 = shift (@fx2)))
-  {
-    my ($f1, $x1) = @$fx1;
-    my ($f2, $x2) = @$fx2;
-
-    die ("Field mismatch $f1 != $f2\n")
-      unless ($f1 eq $f2);
-
-    chomp ($x1); chomp ($x2);
-    if (($x1 !~ m/^\s*$/o) && ($x2 !~ m/^\s*$/o))
-      {
-        my $dx = $x1 - $x2;
-        my $dr = ($x1+$x2 > 0) ? 2*$dx/($x1+$x2) : 0.;
-
-        my $sdx = sprintf ('%17.9e', $dx);
-        my $sdr = sprintf ('%17.9e', $dr);
-
-        $dx = $sdx; $dx = $dx + 0.;
-        $dr = $sdr; $dr = $dr + 0.;
-
-        push @{$x[-1]},
-          sprintf (" $tag1 | %20s | %17.9e  |  %17.9e  |  %17s  |  %17s %s\n", &center ($f1, 20), $x1, $x2, $sdx, $sdr, 
-                   $dr > $opts{'norm-max-diff'} ? '*' : '');
-
-        $nout++ 
-         if ($dr > $opts{'norm-max-diff'});
-
-        if (abs ($dr) > 0)
-          {
-            my $n = &floor ((log (abs ($dr)) / log (10)));
-            $diff{$n}++;
-
-          }
-        else
-          {
-            $zero++;
-          }
-
-        $numb++;
-      }
-    else
-      {
-        push @x, [];
-      }
-  }
-
-printf " $tag1 |                      |%-19s | %-19s | %-19s | %-19s\n", 
-        &center ("NORM(REF)", 19), &center ("NORM(EXP)", 19), 
-        &center ("NORM(REF)-NORM(EXP)", 19), &center ("(NORM(REF)-NORM(EXP))", 19);
-
-printf " $tag1 |                      |%-19s | %-19s | %-19s | %-19s\n", 
-        '', '', '', &center ("/NORM(REF)", 19);
-
-for (my $i = 0; $i <= $#x; $i++)
-  {
-    last unless (@{$x[$i]});
-    print @{$x[$i]};
-  }
-
-print "\n";
-
-
-my $diff_cumul = 0;
-my $perc_cumul = 0;
-for my $n1 (sort { $a <=> $b } keys (%diff))
-  {
-    my $n2 = $n1 + 1;
-    my $diff = $diff{$n1};
-    my $perc = 100 * $diff / $numb;
-    $diff_cumul += $diff;
-    $perc_cumul += $perc;
-    printf (" $tag2 | 1e%+2.2d .. 1e%+2.2d | %3d / %3d | %3d / %3d | %6.2f %%, %6.2f %%\n", $n1, $n2, $diff, $numb, $diff_cumul, $numb, $perc, $perc_cumul);
-  }
-
-if ($nout)
-  {
-    print "\n";
-    my $text = sprintf ("WARNING : SOME NORMS DIFFERENCES ARE OUTSIDE ALLOWED LIMIT OF %6.2f %%\n", 100 * $opts{'norm-max-diff'});
-    print $text x 5;
-  }
-
-print "\n";
-
-my @ot1 = &xobstype ($f1);
-my @ot2 = &xobstype ($f2);
-
-my $not1 = scalar (@ot1);
-my $not2 = scalar (@ot2);
-
-my $not = $not1 > $not2 ? $not1 : $not2;
-
-goto END
-  unless ($not > 0);
-
-my $obs_fmtd  = " %10d |  %13.7e ";
-my $obs_fmtds = " %10s |  %13s ";
-my $obs_fmtp  = " %20d |      %9.4f %% ";
-my $obs_fmtps = " %20s |  %15s ";
-
-printf (" OBS_DIFF | %6s | ", 'Type');
-printf ($obs_fmtds, 'NOBS(REF)', 'JO(REF)');
-print " | ";
-printf ($obs_fmtds, 'NOBS(EXP)', 'JO(EXP)');
-print " | ";
-printf ($obs_fmtps, 'NOBS(EXP)-NOBS(REF)', 'JO(EXP)-JO(EXP)');
-printf ("\n");
-
-my $pot = sub 
-  {
-    my ($ot, $obs_fmt) = @_;
-    $obs_fmt ||= $obs_fmtd;
-    (my $blank = sprintf ($obs_fmt, 0, 0)) =~ s/\S/ /go;
-    if ($ot)
-      {
-        printf ($obs_fmt, $ot->{number}, $ot->{JO});
-      }
-    else
-      {
-        print $blank;
-      }
-  };
-
-my $dot = sub
-  {
-    my ($ot1, $ot2) = @_;
-
-    my $dot = $ot1 && $ot2 ? { number => abs ($ot1->{number} - $ot2->{number}), 
-                               JO => 100 * 2 * abs ($ot1->{JO} - $ot2->{JO}) / abs ($ot1->{JO} + $ot2->{JO}), } : undef;
-    return $dot;
-  };
-
-for my $i (1 .. $not-1)
-  {
-    my $ot1 = $ot1[$i];
-    my $ot2 = $ot2[$i];
-    my $dot12 = $dot->($ot1, $ot2);
-
-    printf (" OBS_DIFF | %6d | ", $i);
-
-    $pot->($ot1);
-    print " | ";
-    $pot->($ot2);
-    print " | ";
-    $pot->($dot12, $obs_fmtp);
-    print "\n";
-
-
-  }
-
-my $jog1 = &xjog ($f1);
-my $jog2 = &xjog ($f2);
-
-if ($jog1 || $jog2)
-  {
-    my $dot12 = $dot->($jog1, $jog2);
-    printf (" OBS_DIFF | GLOBAL | ");
-    $pot->($jog1);
-    print " | ";
-    $pot->($jog2);
-    print " | ";
-    $pot->($dot12, $obs_fmtp);
-    print "\n";
-
-    
-    if ($dot12->{JO} > 100 * $opts{'jo-max-diff'})
-      {
-        my $text = sprintf ("WARNING : GLOBAL JO DIFFERENCE IS OUTSIDE ALLOWED LIMIT OF %12.6f %%\n", 100 * $opts{'jo-max-diff'});
-        print "\n";
-        print $text x 5;
-      }
-
-  }
-
-
-print "\n";
-
-END:
-
-
-
-exit $nout;
-
-