Index: /LMDZ5/trunk/libf/dyn3d_common/disvert.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d_common/disvert.F90	(revision 2044)
+++ /LMDZ5/trunk/libf/dyn3d_common/disvert.F90	(revision 2045)
@@ -43,4 +43,6 @@
   INTEGER l, unit
   REAL dsigmin
+  REAL vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig
+
   REAL alpha, beta, deltaz
   REAL x
@@ -154,5 +156,37 @@
      ap(1)=0.
      ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
-  CASE("strato_custom") 
+  case("strato_correct")
+!==================================================================
+! Fredho 2014/05/18, Saint-Louis du Senegal
+! Cette version de la discretisation strato est corrige au niveau
+! du passage des sig aux ap, bp
+! la version precedente donne un coude dans l'epaisseur des couches
+! vers la tropopause
+!==================================================================
+
+
+     DO l = 1, llm
+        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
+        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
+             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
+     ENDDO
+     dsig = dsig / sum(dsig)
+     sig0(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig0(l) = sig0(l+1) + dsig(l)
+     ENDDO
+     sig=racinesig(sig0)
+
+     bp(1)=1.
+     bp(2:llm)=EXP(1.-1./sig(2: llm)**2)
+     bp(llmp1)=0.
+
+     ap(1)=0.
+     ap(2:llm)=pa*(sig(2:llm)-bp(2:llm))
+     ap(llm+1)=0.
+
+  CASE("strato_custom0") 
+!=======================================================
+! Version Transitoire
     ! custumize strato distribution with specific alpha & beta values and function
     ! depending on llm (experimental and temporary)!
@@ -198,10 +232,73 @@
     ENDIF ! of IF (llm==55.OR.llm==63) ...
     
+
     ! Build sigma distribution
     sig0=EXP(-zz(:)/scaleheight)
     sig0(llm+1)=0.
-    sig=ridders(sig0)
+!    sig=ridders(sig0)
+    sig=racinesig(sig0)
     
     ! Compute ap() and bp()
+    bp(1)=1.
+    bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
+    bp(llm+1)=0.
+    ap=pa*(sig-bp)
+
+  CASE("strato_custom") 
+!===================================================================
+! David Cugnet, François Lott, Lionel Guez, Ehouoarn Millour, Fredho
+! 2014/05
+! custumize strato distribution
+! Al the parameter are given in km assuming a given scalehigh
+    vert_scale=7.     ! scale hight
+    vert_dzmin=0.02   ! width of first layer
+    vert_dzlow=1.     ! dz in the low atmosphere
+    vert_z0low=8.     ! height at which resolution recches dzlow
+    vert_dzmid=3.     ! dz in the mid atmsophere
+    vert_z0mid=70.    ! height at which resolution recches dzmid
+    vert_h_mid=20.    ! width of the transition
+    vert_dzhig=11.    ! dz in the high atmsophere
+    vert_z0hig=80.    ! height at which resolution recches dz
+    vert_h_hig=20.    ! width of the transition
+!===================================================================
+
+    call getin('vert_scale',vert_scale)
+    call getin('vert_dzmin',vert_dzmin)
+    call getin('vert_dzlow',vert_dzlow)
+    call getin('vert_z0low',vert_z0low)
+    CALL getin('vert_dzmid',vert_dzmid)
+    CALL getin('vert_z0mid',vert_z0mid)
+    call getin('vert_h_mid',vert_h_mid)
+    call getin('vert_dzhig',vert_dzhig)
+    call getin('vert_z0hig',vert_z0hig)
+    call getin('vert_h_hig',vert_h_hig)
+
+    scaleheight=vert_scale ! for consistency with further computations
+    ! Build geometrical distribution
+    zz(1)=0.
+    DO l=1,llm
+       z=zz(l)
+       zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*TANH(z/vert_z0low)+                &
+&      (vert_dzmid-vert_dzlow)* &
+&           (TANH((z-vert_z0mid)/vert_h_mid)-TANH((-vert_z0mid)/vert_h_mid)) &
+&     +(vert_dzhig-vert_dzmid-vert_dzlow)*                                  &
+&           (TANH((z-vert_z0hig)/vert_h_hig)-TANH((-vert_z0hig)/vert_h_hig))
+    ENDDO
+
+
+!===================================================================
+! Comment added Fredho 2014/05/18, Saint-Louis, Senegal
+! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H)
+! sig0 is p/p0
+! sig is an intermediate distribution introduce to estimate bp
+! 1.   sig0=exp(-z/H)
+! 2.   inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2)
+! 3.   bp=exp(1-1/sig**2)
+! 4.   ap deduced from  the combination of 2 and 3 so that sig0=ap/p0+bp
+!===================================================================
+
+    sig0=EXP(-zz(:)/vert_scale)
+    sig0(llm+1)=0.
+    sig=racinesig(sig0)
     bp(1)=1.
     bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
@@ -304,4 +401,38 @@
 END FUNCTION ridders
 
+FUNCTION racinesig(sig) RESULT(sg)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Fredho 2014/05/18
+! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
+! Notes:   Uses Newton Raphson search
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: sig(:)
+  REAL              :: sg(SIZE(sig))
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: it, ns, maxit
+  REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
+!-------------------------------------------------------------------------------
+  ns=SIZE(sig); maxit=100
+  c1=Pa/Preff; c2=1.-c1
+  DO l=2,ns-1
+    sg(l)=sig(l)
+    DO it=1,maxit
+       f1=exp(1-1./sg(l)**2)*(1.-c1)
+       sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
+    ENDDO
+!   print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
+  ENDDO
+  sg(1)=1.; sg(ns)=0.
+
+END FUNCTION racinesig
+
+
+
+
 END SUBROUTINE disvert