libf/dyn3d/comvert.f90

module comvert

  use dimens_m, only: llm

  implicit none

  private llm

  real ap(llm+1) ! in Pa
  real bp(llm+1), dpres(llm)
  real presnivs(llm), pa ! in Pa
  real, parameter:: preff = 101325. ! in Pa
  real nivsigs(llm), nivsig(llm+1)

  save

contains

  SUBROUTINE disvert

    ! From dyn3d/disvert.F, v 1.1.1.1 2004/05/19 12:53:05
    ! Auteur : 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.

    use comconst, only: pi

    REAL s(llm+1)
    ! (atmospheric hybrid sigma-pressure coordinate at the interface
    ! between layers "l" and "l-1")

    real ds(llm)
    ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"

    INTEGER l
    REAL alpha, x(llm)

    character(len=7):: s_sampling = "LMD5"
    ! (other allowed values are "param", "strato1" and "strato2")

    real:: h = 7. ! scale height, in km
    ! (used only if "s_sampling" == "param" or "strato1")

    ! These variables are used only in the case 's_sampling == "param"':
    real:: deltaz = 0.04 ! épaisseur de la première 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

    !-----------------------------------------------------------------------

    print *, "Call sequence information: disvert"

    forall (l = 1: llm) nivsigs(l) = REAL(l)
    forall (l = 1: llm + 1) nivsig(l) = REAL(l)

    ! Compute "s":

    s(1) = 1.
    s(llm+1) = 0.

    print *, "Enter namelist 'disvert_nml'."
    read(unit=*, nml=disvert_nml)
    write(unit=*, nml=disvert_nml)

    select case (s_sampling)
    case ("param")
       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
       forall (l = 1: llm) ds(l) &
            = 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2
       ds = ds / sum(ds)

       DO l = llm, 2, -1
          s(l) = s(l+1) + ds(l)
       ENDDO
    case ("strato1")
       ! F. Lott 70 niveaux et plus
       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))
    case ("strato2")
       ! Recommended by F. Lott for a domain including the stratosphere
       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 = ds / sum(ds)

       DO l = llm, 2, -1
          s(l) = s(l+1) + ds(l)
       ENDDO
    case default
       print *, 'Wrong value for "s_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)
       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

  END SUBROUTINE disvert

end module comvert
1	guez	3	module comvert
2
3			use dimens_m, only: llm
4
5			implicit none
6
7			private llm
8
9			real ap(llm+1) ! in Pa
10			real bp(llm+1), dpres(llm)
11			real presnivs(llm), pa ! in Pa
12			real, parameter:: preff = 101325. ! in Pa
13			real nivsigs(llm), nivsig(llm+1)
14
15			save
16
17			contains
18
19			SUBROUTINE disvert
20
21			! From dyn3d/disvert.F, v 1.1.1.1 2004/05/19 12:53:05
22			! Auteur : P. Le Van
23
24			! This procedure sets the vertical grid.
25			! It defines the host variables "ap", "bp", "dpres", "presnivs",
26			! "nivsigs" and "nivsig".
27			! "pa" should be defined before this procedure is called.
28
29			use comconst, only: pi
30
31			REAL s(llm+1)
32			! (atmospheric hybrid sigma-pressure coordinate at the interface
33			! between layers "l" and "l-1")
34
35			real ds(llm)
36			! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"
37
38			INTEGER l
39			REAL alpha, x(llm)
40
41			character(len=7):: s_sampling = "LMD5"
42			! (other allowed values are "param", "strato1" and "strato2")
43
44			real:: h = 7. ! scale height, in km
45			! (used only if "s_sampling" == "param" or "strato1")
46
47			! These variables are used only in the case 's_sampling == "param"':
48			real:: deltaz = 0.04 ! épaisseur de la première couche
49			real:: beta = 1.3 ! facteur d'accroissement en haut
50			real:: k0 = 20. ! nombre de couches dans la transition surface
51			real:: k1 = 1.2 ! nombre de couches dans la transition haute
52
53			REAL ZZ(llm + 1), DZ(llm) ! in km
54
55			namelist /disvert_nml/h, deltaz, beta, k0, k1, s_sampling
56
57			!-----------------------------------------------------------------------
58
59			print *, "Call sequence information: disvert"
60
61			forall (l = 1: llm) nivsigs(l) = REAL(l)
62			forall (l = 1: llm + 1) nivsig(l) = REAL(l)
63
64			! Compute "s":
65
66			s(1) = 1.
67			s(llm+1) = 0.
68
69			print *, "Enter namelist 'disvert_nml'."
70			read(unit=*, nml=disvert_nml)
71			write(unit=*, nml=disvert_nml)
72
73			select case (s_sampling)
74			case ("param")
75			alpha = deltaz / tanh(1./k0) * 2.
76			forall (l = 2: llm) s(l) &
77			= cosh((l - 1) / k0) *(- alpha k0 / h) &
78			* exp(- alpha / h * tanh((llm - k1) / k0) &
79			* beta **(l - 1 - (llm - k1)) / log(beta))
80			case ("LMD5")
81			! Ancienne discrétisation
82			forall (l = 1: llm) ds(l) &
83			= 1. + 7. * SIN(pi * (REAL(l)-0.5) / real(llm+1))**2
84			ds = ds / sum(ds)
85
86			DO l = llm, 2, -1
87			s(l) = s(l+1) + ds(l)
88			ENDDO
89			case ("strato1")
90			! F. Lott 70 niveaux et plus
91			forall (l = 1: llm) dz(l) = 1.56 + TANH(REAL(l - 12) / 5.) &
92			+ TANH(REAL(l - llm) / 10.) / 2.
93
94			zz(1) = 0.
95			DO l = 2, llm + 1
96			zz(l) = zz(l - 1) + dz(l - 1)
97			end DO
98
99			s(2:llm) = (exp(- zz(2:llm) / h) - exp(- zz(llm + 1) / h)) &
100			/ (1. - exp(- zz(llm + 1) / h))
101			case ("strato2")
102			! Recommended by F. Lott for a domain including the stratosphere
103			forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
104
105			ds = (1. + 7. * SIN(x)*2) (1. - tanh(2 * x / pi - 1.))**2 / 4.
106			ds = ds / sum(ds)
107
108			DO l = llm, 2, -1
109			s(l) = s(l+1) + ds(l)
110			ENDDO
111			case default
112			print *, 'Wrong value for "s_sampling"'
113			stop 1
114			END select
115
116			! Calcul de "ap" et "bp" :
117			bp(:llm) = EXP(1. - 1. / s(:llm)**2)
118			bp(llm + 1) = 0.
119			ap = pa * (s - bp)
120
121			forall (l = 1: llm)
122			dpres(l) = bp(l) - bp(l+1)
123			presnivs(l) = 0.5 * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
124			end forall
125
126			END SUBROUTINE disvert
127
128			end module comvert