--- trunk/libf/dyn3d/disvert.f90	2012/09/20 13:00:41	66
+++ trunk/dyn3d/disvert.f	2014/03/12 21:16:36	90
@@ -6,27 +6,27 @@
 
   private llm
 
-  real ap(llm+1), pa ! in Pa
-  real bp(llm+1)
-  real presnivs(llm) ! pressions approximatives des milieux de couches, en Pa
-  real, parameter:: preff = 101325. ! in Pa
-  real nivsigs(llm), nivsig(llm+1)
+  real, save:: ap(llm+1), pa ! in Pa
+  real, save:: bp(llm+1)
+
+  real, save:: presnivs(llm)
+  ! pressions approximatives des milieux de couches, en Pa
 
-  save
+  real, parameter:: preff = 101325. ! in Pa
 
 contains
 
   SUBROUTINE disvert
 
-    ! From dyn3d/disvert.F, v 1.1.1.1 2004/05/19 12:53:05
+    ! From dyn3d/disvert.F, version 1.1.1.1, 2004/05/19 12:53:05
     ! Author: P. Le Van
 
     ! This procedure sets the vertical grid. It defines the host
-    ! variables "ap", "bp", "presnivs", "nivsigs" and "nivsig". "pa"
-    ! should be defined before this procedure is called.
+    ! variables "ap", "bp", "presnivs". "pa" should be defined before
+    ! this procedure is called.
 
     use jumble, only: new_unit
-    use nr_util, only: pi
+    use nr_util, only: pi, assert
     use unit_nml_m, only: unit_nml
 
     ! Local:
@@ -36,41 +36,34 @@
     ! the interface between layers "l" and "l-1"
 
     real ds(llm)
-    ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"
+    ! "ds(l)" : \'epaisseur de la couche "l" dans la coordonn\'ee "s"
 
     INTEGER l, unit
     REAL alpha, x(llm)
 
-    character(len=7):: s_sampling = "tropo"
-    ! (other allowed values are "param", "strato1", "strato" and "read")
+    character(len=7):: vert_sampling = "tropo"
+    ! other allowed values are "param", "strato" and "read"
 
     real:: h = 7. ! scale height, in km
-    ! (used only if "s_sampling" == "param" or "strato1")
+    ! used only if vert_sampling == "param"
 
-    ! These variables are used only in the case 's_sampling == "param"':
-    real:: deltaz = 0.04 ! épaisseur de la première couche
+    ! These variables are used only in the case vert_sampling == "param":
+    real:: deltaz = 0.04 ! \'epaisseur de la premi\`ere couche
     real:: beta = 1.3 ! facteur d'accroissement en haut
     real:: k0 = 20. ! nombre de couches dans la transition surface
     real:: k1 = 1.2 ! nombre de couches dans la transition haute
 
-    REAL ZZ(llm + 1), DZ(llm) ! in km
-
-    namelist /disvert_nml/h, deltaz, beta, k0, k1, s_sampling
+    namelist /disvert_nml/h, deltaz, beta, k0, k1, vert_sampling
 
     !-----------------------------------------------------------------------
 
     print *, "Call sequence information: disvert"
 
-    forall (l = 1: llm) nivsigs(l) = REAL(l)
-    forall (l = 1: llm + 1) nivsig(l) = REAL(l)
-
-    ! Compute "s":
-
     print *, "Enter namelist 'disvert_nml'."
     read(unit=*, nml=disvert_nml)
     write(unit_nml, nml=disvert_nml)
 
-    select case (s_sampling)
+    select case (vert_sampling)
     case ("param")
        s(1) = 1.
        s(llm+1) = 0.
@@ -79,6 +72,7 @@
             = cosh((l - 1) / k0) **(- alpha * k0 / h) &
             * exp(- alpha / h * tanh((llm - k1) / k0) &
             * beta **(l - 1 - (llm - k1)) / log(beta))
+       call compute_ab
     case ("tropo")
        s(1) = 1.
        s(llm+1) = 0.
@@ -89,20 +83,8 @@
        DO l = llm, 2, -1
           s(l) = s(l+1) + ds(l)
        ENDDO
-    case ("strato1")
-       ! F. Lott 70 niveaux et plus
-       s(1) = 1.
-       s(llm+1) = 0.
-       forall (l = 1: llm) dz(l) = 1.56 + TANH(REAL(l - 12) / 5.) &
-               + TANH(REAL(l - llm) / 10.) / 2.
 
-       zz(1) = 0.
-       DO l = 2, llm + 1
-          zz(l) = zz(l - 1) + dz(l - 1)
-       end DO
-
-       s(2:llm) = (exp(- zz(2:llm) / h) - exp(- zz(llm + 1) / h)) &
-            / (1. - exp(- zz(llm + 1) / h))
+       call compute_ab
     case ("strato")
        ! Recommended by F. Lott for a domain including the stratosphere
        s(1) = 1.
@@ -115,33 +97,42 @@
        DO l = llm, 2, -1
           s(l) = s(l+1) + ds(l)
        ENDDO
+
+       call compute_ab
     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.csv", status="old", action="read", &
             position="rewind")
        read(unit, fmt=*) ! skip title line
        do l = 1, llm + 1
-          read(unit, fmt=*) s(l)
+          read(unit, fmt=*) ap(l), bp(l)
        end do
        close(unit)
        ! Quick check:
-       if (s(1) /= 1. .or. s(llm + 1) /= 0. .or. s(1) <= s(2)) then
-          print *, '"s" should be in descending order, from 1 to 0.'
-          stop 1
-       end if
+       call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
+            bp(llm + 1) == 0., "disvert: bad ap or bp values")
     case default
-       print *, 'Wrong value for "s_sampling"'
+       print *, 'Wrong value for "vert_sampling"'
        stop 1
     END select
 
-    ! Calcul de "ap" et "bp" :
-    bp(:llm) = EXP(1. - 1. / s(:llm)**2)
-    bp(llm + 1) = 0.
-    ap = pa * (s - bp)
-
     forall (l = 1: llm) presnivs(l) = 0.5 &
          * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
 
+  contains
+
+    subroutine compute_ab
+
+      ! Calcul de "ap" et "bp" :
+      bp(:llm) = EXP(1. - 1. / s(:llm)**2)
+      bp(llm + 1) = 0.
+      ap = pa * (s - bp)
+
+    end subroutine compute_ab
+
   END SUBROUTINE disvert
 
 end module disvert_m