--- trunk/libf/dyn3d/comvert.f90	2009/07/31 15:18:47	22
+++ trunk/libf/dyn3d/disvert.f90	2012/09/20 13:00:41	66
@@ -1,4 +1,4 @@
-module comvert
+module disvert_m
 
   use dimens_m, only: llm
 
@@ -7,7 +7,7 @@
   private llm
 
   real ap(llm+1), pa ! in Pa
-  real bp(llm+1), dpres(llm)
+  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)
@@ -19,14 +19,17 @@
   SUBROUTINE disvert
 
     ! From dyn3d/disvert.F, v 1.1.1.1 2004/05/19 12:53:05
-    ! Auteur : P. Le Van
+    ! Author: P. Le Van
 
-    ! This procedure sets the vertical grid.
-    ! It defines the host variables "ap", "bp", "dpres", "presnivs",
-    ! "nivsigs" and "nivsig".
-    ! "pa" should be defined before this procedure is called.
+    ! 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.
 
-    use numer_rec, only: pi
+    use jumble, only: new_unit
+    use nr_util, only: pi
+    use unit_nml_m, only: unit_nml
+
+    ! Local:
 
     REAL s(llm+1)
     ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
@@ -35,11 +38,11 @@
     real ds(llm)
     ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"
 
-    INTEGER l
+    INTEGER l, unit
     REAL alpha, x(llm)
 
-    character(len=7):: s_sampling = "LMD5"
-    ! (other allowed values are "param", "strato1" and "strato2")
+    character(len=7):: s_sampling = "tropo"
+    ! (other allowed values are "param", "strato1", "strato" and "read")
 
     real:: h = 7. ! scale height, in km
     ! (used only if "s_sampling" == "param" or "strato1")
@@ -63,24 +66,24 @@
 
     ! Compute "s":
 
-    s(1) = 1.
-    s(llm+1) = 0.
-
     print *, "Enter namelist 'disvert_nml'."
     read(unit=*, nml=disvert_nml)
-    write(unit=*, nml=disvert_nml)
+    write(unit_nml, nml=disvert_nml)
 
     select case (s_sampling)
     case ("param")
+       s(1) = 1.
+       s(llm+1) = 0.
        alpha = deltaz / tanh(1./k0) * 2.
        forall (l = 2: llm) s(l) &
             = cosh((l - 1) / k0) **(- alpha * k0 / h) &
             * exp(- alpha / h * tanh((llm - k1) / k0) &
             * beta **(l - 1 - (llm - k1)) / log(beta))
-    case ("LMD5")
-       ! Ancienne discrétisation
+    case ("tropo")
+       s(1) = 1.
+       s(llm+1) = 0.
        forall (l = 1: llm) ds(l) &
-            = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2
+            = 1. + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2
        ds = ds / sum(ds)
 
        DO l = llm, 2, -1
@@ -88,6 +91,8 @@
        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.
 
@@ -98,16 +103,32 @@
 
        s(2:llm) = (exp(- zz(2:llm) / h) - exp(- zz(llm + 1) / h)) &
             / (1. - exp(- zz(llm + 1) / h))
-    case ("strato2")
+    case ("strato")
        ! Recommended by F. Lott for a domain including the stratosphere
+       s(1) = 1.
+       s(llm+1) = 0.
        forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
 
-       ds = (1. + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.
+       ds = (0.3 + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.
        ds = ds / sum(ds)
 
        DO l = llm, 2, -1
           s(l) = s(l+1) + ds(l)
        ENDDO
+    case("read")
+       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)
+       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
     case default
        print *, 'Wrong value for "s_sampling"'
        stop 1
@@ -118,11 +139,9 @@
     bp(llm + 1) = 0.
     ap = pa * (s - bp)
 
-    forall (l = 1: llm)
-       dpres(l) = bp(l) - bp(l+1)
-       presnivs(l) = 0.5 * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
-    end forall
+    forall (l = 1: llm) presnivs(l) = 0.5 &
+         * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
 
   END SUBROUTINE disvert
 
-end module comvert
+end module disvert_m