Index: trunk/LMDZ.COMMON/libf/dyn3dpar/disvert.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/disvert.F90	(revision 134)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/disvert.F90	(revision 776)
@@ -1,7 +1,11 @@
-! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
+! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
 
 SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
 
   ! Auteur : P. Le Van
+
+  use new_unit_m, only: new_unit
+  use ioipsl, only: getin
+  use assert_m, only: assert
 
   IMPLICIT NONE
@@ -18,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)
@@ -26,20 +31,28 @@
   real zk, zkm1, dzk1, dzk2, k0, k1
 
-  INTEGER l
+  INTEGER l, unit
   REAL dsigmin
   REAL alpha, beta, deltaz
-  INTEGER iostat
   REAL x
   character(len=*),parameter :: modname="disvert"
 
+  character(len=6):: vert_sampling
+  ! (allowed values are "param", "tropo", "strato" and "read")
+
   !-----------------------------------------------------------------------
+
+  print *, "Call sequence information: disvert"
 
   ! default scaleheight is 8km for earth
   scaleheight=8.
 
-  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
+  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
+  call getin('vert_sampling', vert_sampling)
+  print *, 'vert_sampling = ' // vert_sampling
 
-  IF (iostat == 0) THEN
-     ! cas 1 on lit les options dans sigma.def:
+  select case (vert_sampling)
+  case ("param")
+     ! On lit les options dans sigma.def:
+     OPEN(99, file='sigma.def', status='old', form='formatted')
      READ(99, *) scaleheight ! hauteur d'echelle 8.
      READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
@@ -69,31 +82,12 @@
      sig(llm+1)=0.
 
-     DO l = 1, llm
-        dsig(l) = sig(l)-sig(l+1)
-     end DO
-  ELSE
-     if (ok_strato) then
-        if (llm==39) then
-           dsigmin=0.3
-        else if (llm==50) then
-           dsigmin=1.
-        else
-           write(lunout,*) trim(modname), &
-           ' ATTENTION discretisation z a ajuster'
-           dsigmin=1.
-        endif
-        write(lunout,*) trim(modname), &
-        ' Discretisation verticale DSIGMIN=',dsigmin
-     endif
+     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
+     bp(llmp1) = 0.
 
+     ap = pa * (sig - bp)
+  case("tropo")
      DO l = 1, llm
         x = 2*asin(1.) * (l - 0.5) / (llm + 1)
-
-        IF (ok_strato) THEN
-           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
-                *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
-        ELSE
-           dsig(l) = 1.0 + 7.0 * SIN(x)**2
-        ENDIF
+        dsig(l) = 1.0 + 7.0 * SIN(x)**2
      ENDDO
      dsig = dsig / sum(dsig)
@@ -102,5 +96,56 @@
         sig(l) = sig(l+1) + dsig(l)
      ENDDO
-  ENDIF
+
+     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
+        dsigmin=0.3
+     else if (llm==50) then
+        dsigmin=1.
+     else
+        write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
+        dsigmin=1.
+     endif
+     WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
+
+     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)
+     sig(llm+1) = 0.
+     DO l = llm, 1, -1
+        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", &
+          position="rewind")
+     read(unit, fmt=*) ! skip title line
+     do l = 1, llm + 1
+        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)
+  END select
 
   DO l=1, llm
@@ -111,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 '
