libf/dyn3d/disvert.f90

module disvert_m

  use dimens_m, only: llm

  implicit none

  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)

  save

contains

  SUBROUTINE disvert

    ! From dyn3d/disvert.F, v 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.

    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
    ! 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, unit
    REAL alpha, x(llm)

    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")

    ! 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":

    print *, "Enter namelist 'disvert_nml'."
    read(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 ("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
       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
       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))
    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 = (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
    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)

  END SUBROUTINE disvert

end module disvert_m
1	module disvert_m
2
3	use dimens_m, only: llm
4
5	implicit none
6
7	private llm
8
9	real ap(llm+1), pa ! in Pa
10	real bp(llm+1)
11	real presnivs(llm) ! pressions approximatives des milieux de couches, en 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	! Author: P. Le Van
23
24	! This procedure sets the vertical grid. It defines the host
25	! variables "ap", "bp", "presnivs", "nivsigs" and "nivsig". "pa"
26	! should be defined before this procedure is called.
27
28	use jumble, only: new_unit
29	use nr_util, only: pi
30	use unit_nml_m, only: unit_nml
31
32	! Local:
33
34	REAL s(llm+1)
35	! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
36	! the interface between layers "l" and "l-1"
37
38	real ds(llm)
39	! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"
40
41	INTEGER l, unit
42	REAL alpha, x(llm)
43
44	character(len=7):: s_sampling = "tropo"
45	! (other allowed values are "param", "strato1", "strato" and "read")
46
47	real:: h = 7. ! scale height, in km
48	! (used only if "s_sampling" == "param" or "strato1")
49
50	! These variables are used only in the case 's_sampling == "param"':
51	real:: deltaz = 0.04 ! épaisseur de la première couche
52	real:: beta = 1.3 ! facteur d'accroissement en haut
53	real:: k0 = 20. ! nombre de couches dans la transition surface
54	real:: k1 = 1.2 ! nombre de couches dans la transition haute
55
56	REAL ZZ(llm + 1), DZ(llm) ! in km
57
58	namelist /disvert_nml/h, deltaz, beta, k0, k1, s_sampling
59
60	!-----------------------------------------------------------------------
61
62	print *, "Call sequence information: disvert"
63
64	forall (l = 1: llm) nivsigs(l) = REAL(l)
65	forall (l = 1: llm + 1) nivsig(l) = REAL(l)
66
67	! Compute "s":
68
69	print *, "Enter namelist 'disvert_nml'."
70	read(unit=*, nml=disvert_nml)
71	write(unit_nml, nml=disvert_nml)
72
73	select case (s_sampling)
74	case ("param")
75	s(1) = 1.
76	s(llm+1) = 0.
77	alpha = deltaz / tanh(1./k0) * 2.
78	forall (l = 2: llm) s(l) &
79	= cosh((l - 1) / k0) *(- alpha k0 / h) &
80	* exp(- alpha / h * tanh((llm - k1) / k0) &
81	* beta **(l - 1 - (llm - k1)) / log(beta))
82	case ("tropo")
83	s(1) = 1.
84	s(llm+1) = 0.
85	forall (l = 1: llm) ds(l) &
86	= 1. + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2
87	ds = ds / sum(ds)
88
89	DO l = llm, 2, -1
90	s(l) = s(l+1) + ds(l)
91	ENDDO
92	case ("strato1")
93	! F. Lott 70 niveaux et plus
94	s(1) = 1.
95	s(llm+1) = 0.
96	forall (l = 1: llm) dz(l) = 1.56 + TANH(REAL(l - 12) / 5.) &
97	+ TANH(REAL(l - llm) / 10.) / 2.
98
99	zz(1) = 0.
100	DO l = 2, llm + 1
101	zz(l) = zz(l - 1) + dz(l - 1)
102	end DO
103
104	s(2:llm) = (exp(- zz(2:llm) / h) - exp(- zz(llm + 1) / h)) &
105	/ (1. - exp(- zz(llm + 1) / h))
106	case ("strato")
107	! Recommended by F. Lott for a domain including the stratosphere
108	s(1) = 1.
109	s(llm+1) = 0.
110	forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
111
112	ds = (0.3 + 7. * SIN(x)*2) (1. - tanh(2 * x / pi - 1.))**2 / 4.
113	ds = ds / sum(ds)
114
115	DO l = llm, 2, -1
116	s(l) = s(l+1) + ds(l)
117	ENDDO
118	case("read")
119	call new_unit(unit)
120	open(unit, file="hybrid.csv", status="old", action="read", &
121	position="rewind")
122	read(unit, fmt=*) ! skip title line
123	do l = 1, llm + 1
124	read(unit, fmt=*) s(l)
125	end do
126	close(unit)
127	! Quick check:
128	if (s(1) /= 1. .or. s(llm + 1) /= 0. .or. s(1) <= s(2)) then
129	print *, '"s" should be in descending order, from 1 to 0.'
130	stop 1
131	end if
132	case default
133	print *, 'Wrong value for "s_sampling"'
134	stop 1
135	END select
136
137	! Calcul de "ap" et "bp" :
138	bp(:llm) = EXP(1. - 1. / s(:llm)**2)
139	bp(llm + 1) = 0.
140	ap = pa * (s - bp)
141
142	forall (l = 1: llm) presnivs(l) = 0.5 &
143	* (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
144
145	END SUBROUTINE disvert
146
147	end module disvert_m