Index: /LMDZ5/trunk/libf/dyn3d/disvert.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/disvert.F90	(revision 1644)
+++ /LMDZ5/trunk/libf/dyn3d/disvert.F90	(revision 1645)
@@ -7,4 +7,5 @@
   use new_unit_m, only: new_unit
   use ioipsl, only: getin
+  use assert_m, only: assert
 
   IMPLICIT NONE
@@ -21,5 +22,6 @@
 
   real,intent(in) :: pa, preff
-  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: ap(llmp1) ! in Pa
+  real, intent(out):: bp(llmp1)
   real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
   real,intent(out) :: presnivs(llm)
@@ -79,4 +81,9 @@
 
      sig(llm+1)=0.
+
+     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
+     bp(llmp1) = 0.
+
+     ap = pa * (sig - bp)
   case("tropo")
      DO l = 1, llm
@@ -89,4 +96,11 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
   case("strato")
      if (llm==39) then
@@ -110,5 +124,15 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
   case("read")
+     ! Read "ap" and "bp". First line is skipped (title line). "ap"
+     ! should be in Pa. First couple of values should correspond to
+     ! the surface, that is : "bp" should be in descending order.
      call new_unit(unit)
      open(unit, file="hybrid.txt", status="old", action="read", &
@@ -116,7 +140,9 @@
      read(unit, fmt=*) ! skip title line
      do l = 1, llm + 1
-        read(unit, fmt=*) sig(l)
+        read(unit, fmt=*) ap(l), bp(l)
      end do
      close(unit)
+     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
+          bp(llm + 1) == 0., "disvert: bad ap or bp values")
   case default
      call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
@@ -130,19 +156,4 @@
      nivsig(l)= REAL(l)
   ENDDO
-
-  ! .... Calculs de ap(l) et de bp(l) ....
-  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
-
-  bp(llmp1) = 0.
-
-  DO l = 1, llm
-     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
-     ap(l) = pa * ( sig(l) - bp(l) )
-  ENDDO
-
-  bp(1)=1.
-  ap(1)=0.
-
-  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
   write(lunout, *)  trim(modname),': BP '
Index: /LMDZ5/trunk/libf/dyn3dpar/disvert.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3dpar/disvert.F90	(revision 1644)
+++ /LMDZ5/trunk/libf/dyn3dpar/disvert.F90	(revision 1645)
@@ -7,4 +7,5 @@
   use new_unit_m, only: new_unit
   use ioipsl, only: getin
+  use assert_m, only: assert
 
   IMPLICIT NONE
@@ -21,5 +22,6 @@
 
   real,intent(in) :: pa, preff
-  real,intent(out) :: ap(llmp1), bp(llmp1)
+  real,intent(out) :: ap(llmp1) ! in Pa
+  real, intent(out):: bp(llmp1)
   real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
   real,intent(out) :: presnivs(llm)
@@ -79,4 +81,9 @@
 
      sig(llm+1)=0.
+
+     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
+     bp(llmp1) = 0.
+
+     ap = pa * (sig - bp)
   case("tropo")
      DO l = 1, llm
@@ -89,4 +96,11 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
   case("strato")
      if (llm==39) then
@@ -110,5 +124,15 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
+
+     bp(1)=1.
+     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
+     bp(llmp1) = 0.
+
+     ap(1)=0.
+     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
   case("read")
+     ! Read "ap" and "bp". First line is skipped (title line). "ap"
+     ! should be in Pa. First couple of values should correspond to
+     ! the surface, that is : "bp" should be in descending order.
      call new_unit(unit)
      open(unit, file="hybrid.txt", status="old", action="read", &
@@ -116,7 +140,9 @@
      read(unit, fmt=*) ! skip title line
      do l = 1, llm + 1
-        read(unit, fmt=*) sig(l)
+        read(unit, fmt=*) ap(l), bp(l)
      end do
      close(unit)
+     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
+          bp(llm + 1) == 0., "disvert: bad ap or bp values")
   case default
      call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
@@ -130,19 +156,4 @@
      nivsig(l)= REAL(l)
   ENDDO
-
-  ! .... Calculs de ap(l) et de bp(l) ....
-  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
-
-  bp(llmp1) = 0.
-
-  DO l = 1, llm
-     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
-     ap(l) = pa * ( sig(l) - bp(l) )
-  ENDDO
-
-  bp(1)=1.
-  ap(1)=0.
-
-  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
 
   write(lunout, *)  trim(modname),': BP '