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