Index: trunk/UTIL/z2sig/README
===================================================================
--- trunk/UTIL/z2sig/README	(revision 3651)
+++ trunk/UTIL/z2sig/README	(revision 3651)
@@ -0,0 +1,3 @@
+Scripts to build a new z2sig.def using an input z2sig.def, but with a different number of layers.
+
+z2sig.def_47 is an example for Pluto.
Index: trunk/UTIL/z2sig/build_sig.py
===================================================================
--- trunk/UTIL/z2sig/build_sig.py	(revision 3651)
+++ trunk/UTIL/z2sig/build_sig.py	(revision 3651)
@@ -0,0 +1,193 @@
+#!/usr/bin/python3
+from	numpy		      import	*
+import  numpy                 as        np
+import  matplotlib.pyplot     as        mpl
+from matplotlib.cm import get_cmap
+import pylab
+from matplotlib import ticker
+import matplotlib.colors as colors
+import datetime
+from matplotlib import pyplot
+
+######################################################################
+# Fonction  read  file
+######################################################################
+def readfile(name):
+# test reading lire texte txt
+    myfile = open(name, 'r')
+    mylines=myfile.readlines()
+    nbline=size(mylines)
+    nbcolumn=len(mylines[0].split())
+    print('nbline, nbcol:',nbline, nbcolumn)
+    data=np.zeros(nbline-1,dtype='f')
+    i=0
+    for line in mylines:
+        s=str.split(line)
+        if i==0:
+          hh=float(s[0])
+        else:
+          data[i-1]=float(s[0])
+        i=i+1
+    return data,hh
+
+######################################################################
+# Input Parameters
+######################################################################
+nbl=47
+data,hh=readfile('z2sig.def_47')
+# nb of levels
+clev=np.linspace(1,nbl,nbl)
+ps=1.1
+# New levels
+mynbl=201
+n=float(mynbl)
+plot = True
+plot = False
+
+#####################
+# Get sigma levels
+#####################
+sig=np.ones(nbl,dtype='f')
+for i in range(nbl-1):
+ sig[i+1]=0.5*(exp(-data[i+1]/hh)+exp(-data[i]/hh))
+
+#####################
+# Get non-hybrid pressure levels
+#####################
+ap=0
+bp=sig
+pp=ap+bp*1.1
+
+#####################
+# Get hybrid pressure levels
+#####################
+def sig_hybrid(sig,pa,preff):
+#c     Subroutine utilisee pour calculer des valeurs de sigma modifie
+#c     pour conserver les coordonnees verticales decrites dans
+#c     esasig.def/z2sig.def lors du passage en coordonnees hybrides
+#c     F. Forget 2002
+#c     Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
+#c     L'objectif est de calculer newsig telle que
+#c       (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
+#c     Cela ne se resoud pas analytiquement:
+#c     => on resoud par iterration bourrine
+#c     ----------------------------------------------
+#c     Information  : where exp(1-1./x**2) become << x
+#c           x      exp(1-1./x**2) /x
+#c           1           1
+#c           0.68       0.5
+#c           0.391      1.E-2
+#c        => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25
+
+  newsig = sig
+  x1=0
+  x2=1
+  if (sig>=1):
+     newsig = sig
+  elif (sig*preff/pa>=0.25):
+     for j in range(9999):  # nombre d'iteration max
+       F=((1-pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
+       if (F>1):
+          x2 = newsig
+          newsig=(x1+newsig)*0.5
+       else:
+          x1=newsig
+          newsig=(x2+newsig)*0.5
+       #Test : on arete lorsque on approxime sig a moins de 0.01 m pres
+       #(en pseudo altiude) :
+       #if (abs(10.*log(F))<1.e-5):
+       # break
+
+  else:  # if (sig*preff/pa.le.0.25)
+      newsig= sig*preff/pa
+
+  return newsig
+#----------------------------------
+
+# hybrid pressure:
+nsig=np.zeros(nbl,dtype='f')
+pa=0.5 #ps
+preff=2 #ps
+for i in range(nbl):
+  nsig[i]=sig_hybrid(sig[i],pa,preff)
+
+bps=exp(1-1/(nsig)**2)
+aps=pa*(nsig-bps)
+
+pps=aps+bps*ps
+
+#####################
+# Get Pseudo altitude
+#####################
+# non hybrid
+ph=-hh*log(pp/ps)
+dzph=-(ph[0:-1]-ph[1:])
+# hybrid
+phs=-hh*log(pps/ps)
+dzphs=-(phs[0:-1]-phs[1:])
+
+#####################
+# New choice delta-zh
+#####################
+myclev=np.linspace(1,mynbl,mynbl)
+xx=myclev[0:mynbl-1] # xaxis for layer thickness
+xx1=myclev[0:mynbl]
+# first and last level difference in m
+l0=0.0025
+ln=4.6
+# Coefficient for the exp function
+c=0.11
+a=1/n*(np.log(ln/c)-np.log(l0/c))
+b=-1/a*np.log(l0/c)
+# Dz difference altitude between 2 layers :
+mydz=exp(a*(xx-b))*c
+
+# New pseudo altitude myph
+myph=np.zeros(mynbl,dtype='f')
+myph[0]=mydz[0]
+print(myph[0])
+for i in range(mynbl-1):
+ myph[i+1]=mydz[i]+myph[i]
+ print(myph[i+1])
+
+to_file = np.concatenate([[hh], myph])
+output = f"z2sig.def_{mynbl}"
+np.savetxt(output,to_file, fmt="%g")
+print(f"Saved to {output}")
+
+# New Ps levels
+mypp=ps*exp(-myph/hh)
+
+
+if plot:
+    ### Plot pressure levels
+    mpl.figure(figsize=(50/2.54, 30/2.54),facecolor='w')
+    ax=mpl.subplot(131)
+    mpl.plot(clev,pp,'k-',marker='*',label='P ini')
+    mpl.plot(clev,pps,'r-',marker='*',label='P ini hybrid')
+    mpl.plot(myclev,mypp,'b-',marker='*',label='P new')
+    mpl.gca().invert_yaxis()
+    pyplot.yscale('log')
+    mpl.grid()
+
+    ### Plot difference altitude between 2 layers
+    ax=mpl.subplot(132)
+    mpl.plot(clev[0:nbl-1],dzph,'k-',marker='*',label='dz ini')
+    mpl.plot(clev[0:nbl-1],dzphs,'r-',marker='*',label='dz ini hybrid')
+    mpl.plot(xx,mydz,'b-',marker='*',label='dz new')
+    pyplot.yscale('log')
+    mpl.grid()
+
+    ### pseudo altitude
+    ax=mpl.subplot(133)
+    mpl.plot(clev[0:nbl],ph,'k-',marker='*',label='sig ini')
+    mpl.plot(clev[0:nbl],phs,'r-',marker='*',label='sig ini hybrid')
+    mpl.plot(xx1,myph,'b-',marker='*',label='sig new')
+    #pyplot.yscale('log')
+    mpl.grid()
+
+    mpl.show()
+
+
+
+
Index: trunk/UTIL/z2sig/build_sig_uneven.py
===================================================================
--- trunk/UTIL/z2sig/build_sig_uneven.py	(revision 3651)
+++ trunk/UTIL/z2sig/build_sig_uneven.py	(revision 3651)
@@ -0,0 +1,27 @@
+#!/opt/canopy-1.3.0/Canopy_64bit/User/bin/python
+import numpy as np
+
+def compute_uneven_sigma(num_levels, N_scale_heights, surf_res, exponent, zero_top ):
+    """
+    Construct an initial array of sigma based on the number of levels, an exponent
+    Args:
+        num_levels: the number of levels
+        N_scale_heights: the number of scale heights to the top of the model (e.g scale_heights =12.5 ~102km assuming 8km scale height)
+        surf_res: the resolution at the surface
+        exponent: an exponent to increase th thickness of the levels
+        zero_top: if True, force the top pressure boundary (in N=0) to 0 Pa
+    Returns:
+        b: an array of sigma layers
+
+    """
+    b=np.zeros(int(num_levels)+1)
+    for k in range(0,int(num_levels)):
+        zeta = 1.-k/np.float_(num_levels) #zeta decreases with k
+        z  = surf_res*zeta + (1.0 - surf_res)*(zeta**exponent)
+        b[k] = np.exp(-z*N_scale_heights)
+    b[-1] = 1.0
+    if(zero_top):  b[0] = 0.0
+    return b
+
+sig=compute_uneven_sigma(47, 10, 5, 1, 0 )
+print(sig)
Index: trunk/UTIL/z2sig/z2sig.def_47
===================================================================
--- trunk/UTIL/z2sig/z2sig.def_47	(revision 3651)
+++ trunk/UTIL/z2sig/z2sig.def_47	(revision 3651)
@@ -0,0 +1,48 @@
+18.00000
+0.010
+0.015
+0.025
+0.040
+0.060
+0.090
+0.1200
+0.1500
+0.1800
+0.2200
+0.2800
+0.3400
+0.5000
+0.70000
+0.90000
+1.20000
+1.50000
+1.80000
+2.10000
+2.50000
+3.00000
+3.60000
+5.00000
+6.50000
+8.00000
+9.50000
+11.0000
+13.0000
+16.0000
+19.0000
+22.0000
+25.0000
+29.0000
+32.0000
+36.0000
+40.0000
+45.0000
+50.0000
+55.0000
+60.0000
+66.0000
+72.0000
+80.0000
+90.0000
+105.000
+120.000
+140.000