1 |
module disvert_m |
module disvert_m |
2 |
|
|
3 |
use dimens_m, only: llm |
use dimensions, only: llm |
4 |
|
|
5 |
implicit none |
implicit none |
6 |
|
|
7 |
private llm |
private llm, hybrid, funcd, y, ya, compute_ab |
8 |
|
|
9 |
real ap(llm+1), pa ! in Pa |
real, save:: ap(llm+1), pa ! in Pa |
10 |
real bp(llm+1) |
real, save:: bp(llm+1) |
11 |
real presnivs(llm) ! pressions approximatives des milieux de couches, en Pa |
|
12 |
real, parameter:: preff = 101325. ! in Pa |
REAL s(llm+1) |
13 |
real nivsigs(llm), nivsig(llm+1) |
! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at |
14 |
|
! half-level, between layers "l" and "l-1" |
15 |
|
|
16 |
save |
real, save:: presnivs(llm) |
17 |
|
! approximate full level pressure for a reference surface pressure, in Pa |
18 |
|
|
19 |
|
real, parameter:: preff = 101325. ! in Pa |
20 |
|
real y, ya ! for the hybrid function |
21 |
|
|
22 |
contains |
contains |
23 |
|
|
24 |
SUBROUTINE disvert |
SUBROUTINE disvert |
25 |
|
|
26 |
! 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 |
27 |
! Author: P. Le Van |
! Author: P. Le Van |
28 |
|
|
29 |
! This procedure sets the vertical grid. It defines the host |
! This procedure sets the vertical grid. It defines the host |
30 |
! variables "ap", "bp", "presnivs", "nivsigs" and "nivsig". "pa" |
! variables "ap", "bp", "presnivs". "pa" should be defined before |
31 |
! should be defined before this procedure is called. |
! this procedure is called. |
32 |
|
|
33 |
use jumble, only: new_unit |
use jumble, only: read_column, new_unit |
34 |
use nr_util, only: pi |
use nr_util, only: pi, assert |
35 |
use unit_nml_m, only: unit_nml |
use unit_nml_m, only: unit_nml |
36 |
|
|
37 |
! Local: |
! Local: |
38 |
|
|
|
REAL s(llm+1) |
|
|
! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at |
|
|
! the interface between layers "l" and "l-1" |
|
|
|
|
39 |
real ds(llm) |
real ds(llm) |
40 |
! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s" |
! "ds(l)" : \'epaisseur de la couche "l" dans la coordonn\'ee "s" |
41 |
|
|
42 |
INTEGER l, unit |
INTEGER l, unit |
43 |
REAL alpha, x(llm) |
REAL x(llm) |
44 |
|
real:: dsigmin = 1. |
45 |
character(len=7):: s_sampling = "tropo" |
real zz(llm) ! in km |
46 |
! (other allowed values are "param", "strato1", "strato" and "read") |
|
47 |
|
character(len=20):: vert_sampling = "tropo" |
48 |
real:: h = 7. ! scale height, in km |
! Allowed values: "tropo", "strato_custom", "strato", |
49 |
! (used only if "s_sampling" == "param" or "strato1") |
! "read_hybrid", "read_pressure". |
50 |
|
|
51 |
! These variables are used only in the case 's_sampling == "param"': |
! These variables are used only in the case vert_sampling == |
52 |
real:: deltaz = 0.04 ! épaisseur de la première couche |
! "strato_custom", and all are in km: |
53 |
real:: beta = 1.3 ! facteur d'accroissement en haut |
real:: vert_scale = 7. ! scale height |
54 |
real:: k0 = 20. ! nombre de couches dans la transition surface |
real:: vert_dzmin = 0.017 ! width of first layer |
55 |
real:: k1 = 1.2 ! nombre de couches dans la transition haute |
real:: vert_dzlow = 1. ! dz in the low atmosphere |
56 |
|
real:: vert_z0low = 8.7 ! height at which resolution reaches dzlow |
57 |
REAL ZZ(llm + 1), DZ(llm) ! in km |
real:: vert_dzmid = 2. ! dz in the mid atmosphere |
58 |
|
real:: vert_z0mid = 70. ! height at which resolution reaches dzmid |
59 |
namelist /disvert_nml/h, deltaz, beta, k0, k1, s_sampling |
real:: vert_h_mid = 20. ! width of the transition |
60 |
|
real:: vert_dzhig = 11. ! dz in the high atmosphere |
61 |
|
real:: vert_z0hig = 75. ! height at which resolution reaches dz |
62 |
|
real:: vert_h_hig = 20. ! width of the transition |
63 |
|
|
64 |
|
real, allocatable:: p(:) ! (2:llm or llm + 1) pressure (in hPa) |
65 |
|
|
66 |
|
namelist /disvert_nml/vert_sampling, vert_scale, vert_dzmin, vert_dzlow, & |
67 |
|
vert_z0low, vert_dzmid, vert_z0mid, vert_h_mid, vert_dzhig, & |
68 |
|
vert_z0hig, vert_h_hig, dsigmin |
69 |
|
|
70 |
!----------------------------------------------------------------------- |
!----------------------------------------------------------------------- |
71 |
|
|
72 |
print *, "Call sequence information: disvert" |
print *, "Call sequence information: disvert" |
73 |
|
|
74 |
forall (l = 1: llm) nivsigs(l) = REAL(l) |
write(unit=*, nml=disvert_nml) |
|
forall (l = 1: llm + 1) nivsig(l) = REAL(l) |
|
|
|
|
|
! Compute "s": |
|
|
|
|
75 |
print *, "Enter namelist 'disvert_nml'." |
print *, "Enter namelist 'disvert_nml'." |
76 |
read(unit=*, nml=disvert_nml) |
read(unit=*, nml=disvert_nml) |
77 |
write(unit_nml, nml=disvert_nml) |
write(unit_nml, nml=disvert_nml) |
78 |
|
|
79 |
select case (s_sampling) |
s(1) = 1. |
80 |
case ("param") |
s(llm+1) = 0. |
81 |
s(1) = 1. |
|
82 |
s(llm+1) = 0. |
select case (vert_sampling) |
83 |
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)) |
|
84 |
case ("tropo") |
case ("tropo") |
85 |
s(1) = 1. |
! with llm = 19 and dsigmin = 1 for CMIP 3 |
86 |
s(llm+1) = 0. |
|
87 |
forall (l = 1: llm) ds(l) & |
forall (l = 1: llm) ds(l) & |
88 |
= 1. + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2 |
= dsigmin + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2 |
89 |
ds = ds / sum(ds) |
ds = ds / sum(ds) |
90 |
|
|
91 |
DO l = llm, 2, -1 |
DO l = llm, 2, -1 |
92 |
s(l) = s(l+1) + ds(l) |
s(l) = s(l+1) + ds(l) |
93 |
ENDDO |
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. |
|
94 |
|
|
95 |
zz(1) = 0. |
call compute_ab |
|
DO l = 2, llm + 1 |
|
|
zz(l) = zz(l - 1) + dz(l - 1) |
|
|
end DO |
|
96 |
|
|
|
s(2:llm) = (exp(- zz(2:llm) / h) - exp(- zz(llm + 1) / h)) & |
|
|
/ (1. - exp(- zz(llm + 1) / h)) |
|
97 |
case ("strato") |
case ("strato") |
98 |
! Recommended by F. Lott for a domain including the stratosphere |
! with llm = 39 and dsigmin = 0.3 for CMIP5 |
99 |
s(1) = 1. |
|
|
s(llm+1) = 0. |
|
100 |
forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1) |
forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1) |
101 |
|
|
102 |
ds = (0.3 + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4. |
ds = (dsigmin + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4. |
103 |
ds = ds / sum(ds) |
ds = ds / sum(ds) |
104 |
|
|
105 |
DO l = llm, 2, -1 |
DO l = llm, 2, -1 |
106 |
s(l) = s(l+1) + ds(l) |
s(l) = s(l+1) + ds(l) |
107 |
ENDDO |
ENDDO |
108 |
case("read") |
|
109 |
|
call compute_ab |
110 |
|
|
111 |
|
case ("strato_custom") |
112 |
|
! with llm = 79 for CMIP 6 |
113 |
|
|
114 |
|
zz(1) = 0. |
115 |
|
DO l = 1, llm - 1 |
116 |
|
zz(l + 1) = zz(l) + vert_dzmin + vert_dzlow & |
117 |
|
* TANH(zz(l) / vert_z0low) + (vert_dzmid - vert_dzlow) & |
118 |
|
* (TANH((zz(l) - vert_z0mid) / vert_h_mid) & |
119 |
|
- TANH(- vert_z0mid / vert_h_mid)) & |
120 |
|
+ (vert_dzhig - vert_dzmid - vert_dzlow) & |
121 |
|
* (TANH((zz(l) - vert_z0hig) / vert_h_hig) & |
122 |
|
- TANH(- vert_z0hig / vert_h_hig)) |
123 |
|
ENDDO |
124 |
|
|
125 |
|
allocate(p(2: llm)) |
126 |
|
p = preff * EXP(- zz(2:) / vert_scale) |
127 |
|
ya = pa / preff |
128 |
|
s(2: llm) = hybrid(p) |
129 |
|
|
130 |
|
call compute_ab |
131 |
|
|
132 |
|
case("read_hybrid") |
133 |
|
! Read "ap" and "bp". First line is skipped (title line). "ap" |
134 |
|
! should be in Pa. First couple of values should correspond to |
135 |
|
! the surface, that is : "bp" should be in descending order. |
136 |
call new_unit(unit) |
call new_unit(unit) |
137 |
open(unit, file="hybrid.csv", status="old", action="read", & |
open(unit, file="hybrid.csv", status="old", action="read", & |
138 |
position="rewind") |
position="rewind") |
139 |
read(unit, fmt=*) ! skip title line |
read(unit, fmt=*) ! skip title line |
140 |
do l = 1, llm + 1 |
do l = 1, llm + 1 |
141 |
read(unit, fmt=*) s(l) |
read(unit, fmt=*) ap(l), bp(l) |
142 |
end do |
end do |
143 |
close(unit) |
close(unit) |
144 |
! Quick check: |
! Quick check: |
145 |
if (s(1) /= 1. .or. s(llm + 1) /= 0. .or. s(1) <= s(2)) then |
call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., & |
146 |
print *, '"s" should be in descending order, from 1 to 0.' |
bp(llm + 1) == 0., "disvert: bad ap or bp values") |
147 |
stop 1 |
s(2: llm) = ap(2: llm) / pa + bp(2: llm) |
148 |
end if |
|
149 |
|
case("read_pressure") |
150 |
|
! Read pressure values, in Pa, in descending order, from preff |
151 |
|
! to 0. First line is skipped (title line). |
152 |
|
call read_column("pressure.txt", p, first=2) |
153 |
|
call assert(size(p) == llm + 1, "disvert: bad number of pressure values") |
154 |
|
! Quick check: |
155 |
|
call assert(p(1) == preff, p(llm + 1) == 0., & |
156 |
|
"disvert: bad pressure values") |
157 |
|
ya = pa / preff |
158 |
|
s(2: llm) = hybrid(p(2: llm)) |
159 |
|
call compute_ab |
160 |
|
|
161 |
case default |
case default |
162 |
print *, 'Wrong value for "s_sampling"' |
print *, 'Wrong value for "vert_sampling"' |
163 |
stop 1 |
stop 1 |
|
END select |
|
164 |
|
|
165 |
! Calcul de "ap" et "bp" : |
END select |
|
bp(:llm) = EXP(1. - 1. / s(:llm)**2) |
|
|
bp(llm + 1) = 0. |
|
|
ap = pa * (s - bp) |
|
166 |
|
|
167 |
forall (l = 1: llm) presnivs(l) = 0.5 & |
forall (l = 1: llm) presnivs(l) = 0.5 & |
168 |
* (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff) |
* (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff) |
169 |
|
|
170 |
END SUBROUTINE disvert |
END SUBROUTINE disvert |
171 |
|
|
172 |
|
!********************************************************** |
173 |
|
|
174 |
|
subroutine compute_ab |
175 |
|
|
176 |
|
! Calcul de "ap" et "bp". |
177 |
|
|
178 |
|
where (s >= 1. / sqrt(1. - log(tiny(0.)))) |
179 |
|
bp = exp(1. - 1. / s**2) |
180 |
|
elsewhere |
181 |
|
bp = 0. |
182 |
|
end where |
183 |
|
|
184 |
|
ap = pa * (s - bp) |
185 |
|
|
186 |
|
end subroutine compute_ab |
187 |
|
|
188 |
|
!********************************************************** |
189 |
|
|
190 |
|
function hybrid(p) |
191 |
|
|
192 |
|
! This procedure computes the hybrid sigma-pressure coordinate |
193 |
|
! from pressure values, assuming some reference surface |
194 |
|
! pressure. The procedure assumes, and does not check, that |
195 |
|
! pressure values are given in descending order. |
196 |
|
|
197 |
|
use numer_rec_95, only: rtsafe |
198 |
|
|
199 |
|
real, intent(in):: p(:) ! pressure (in hPa) |
200 |
|
real hybrid(size(p)) ! hybrid sigma-pressure coordinate |
201 |
|
|
202 |
|
! Local: |
203 |
|
integer l |
204 |
|
|
205 |
|
!------------------------------------------------------- |
206 |
|
|
207 |
|
y = p(1) / preff |
208 |
|
hybrid(1) = rtsafe(funcd, x1 = 0., x2 = 1., xacc = 1e-4) |
209 |
|
|
210 |
|
do l = 2, size(p) |
211 |
|
y = p(l) / preff |
212 |
|
! Assuming descending order in pressure: |
213 |
|
hybrid(l) = rtsafe(funcd, x1 = 0., x2 = hybrid(l - 1), & |
214 |
|
xacc = hybrid(l - 1) * 1e-4) |
215 |
|
end do |
216 |
|
|
217 |
|
end function hybrid |
218 |
|
|
219 |
|
!********************************************************** |
220 |
|
|
221 |
|
SUBROUTINE funcd(s, fval, fderiv) |
222 |
|
|
223 |
|
REAL, INTENT(IN):: s |
224 |
|
REAL, INTENT(OUT):: fval, fderiv |
225 |
|
|
226 |
|
! Local: |
227 |
|
real b |
228 |
|
|
229 |
|
!------------------------------------ |
230 |
|
|
231 |
|
if (s**3 > 1. / huge(0.)) then |
232 |
|
b = exp(1. - 1. / s**2) |
233 |
|
fval = ya * s + (1. - ya) * b - y |
234 |
|
fderiv = ya + 2 * (1. - ya) * b / s**3 |
235 |
|
else |
236 |
|
fval = ya * s - y |
237 |
|
fderiv = ya |
238 |
|
end if |
239 |
|
|
240 |
|
END SUBROUTINE funcd |
241 |
|
|
242 |
end module disvert_m |
end module disvert_m |