GCC Code Coverage Report

Directory: ./
File: dyn3d_common/disvert.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 46 170 27.1%
Branches: 21 95 22.1%
Line Branch Exec Source
1 ! $Id: disvert.F90 2786 2017-01-30 10:24:47Z emillour $
2
3 52 SUBROUTINE disvert()
4
5 use ioipsl, only: getin
6 use new_unit_m, only: new_unit
7 use assert_m, only: assert
8 USE comvert_mod, ONLY: ap, bp, aps, bps, nivsigs, nivsig, dpres, presnivs, &
9 pseudoalt, pa, preff, scaleheight
10 USE logic_mod, ONLY: ok_strato
11
12 IMPLICIT NONE
13
14 include "dimensions.h"
15 include "paramet.h"
16 include "iniprint.h"
17
18 !-------------------------------------------------------------------------------
19 ! Purpose: Vertical distribution functions for LMDZ.
20 ! Triggered by the levels number llm.
21 !-------------------------------------------------------------------------------
22 ! Read in "comvert_mod":
23
24 ! pa !--- vertical coordinate is close to a PRESSURE COORDINATE FOR P
25 ! < 0.3 * pa (relative variation of p on a model level is < 0.1 %)
26
27 ! preff !--- REFERENCE PRESSURE (101325 Pa)
28 ! Written in "comvert_mod":
29 ! ap(llm+1), bp(llm+1) !--- Ap, Bp HYBRID COEFFICIENTS AT INTERFACES
30 ! aps(llm), bps(llm) !--- Ap, Bp HYBRID COEFFICIENTS AT MID-LAYERS
31 ! dpres(llm) !--- PRESSURE DIFFERENCE FOR EACH LAYER
32 ! presnivs(llm) !--- PRESSURE AT EACH MID-LAYER
33 ! scaleheight !--- VERTICAL SCALE HEIGHT (Earth: 8kms)
34 ! nivsig(llm+1) !--- SIGMA INDEX OF EACH LAYER INTERFACE
35 ! nivsigs(llm) !--- SIGMA INDEX OF EACH MID-LAYER
36 !-------------------------------------------------------------------------------
37 ! Local variables:
38 REAL sig(llm+1), dsig(llm)
39 REAL sig0(llm+1), zz(llm+1)
40 REAL zk, zkm1, dzk1, dzk2, z, k0, k1
41
42 INTEGER l, unit
43 REAL dsigmin
44 REAL vert_scale,vert_dzmin,vert_dzlow,vert_z0low,vert_dzmid,vert_z0mid,vert_h_mid,vert_dzhig,vert_z0hig,vert_h_hig
45
46 REAL alpha, beta, deltaz
47 REAL x
48 character(len=*),parameter :: modname="disvert"
49
50 character(len=24):: vert_sampling
51 ! (allowed values are "param", "tropo", "strato" and "read")
52
53 !-----------------------------------------------------------------------
54
55 1 WRITE(lunout,*) TRIM(modname)//" starts"
56
57 ! default scaleheight is 8km for earth
58 1 scaleheight=8.
59
60
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
61 1 call getin('vert_sampling', vert_sampling)
62 1 WRITE(lunout,*) TRIM(modname)//' vert_sampling = ' // vert_sampling
63
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (llm==39 .and. vert_sampling=="strato") then
64 1 dsigmin=0.3 ! Vieille option par défaut pour CMIP5
65 else
66 ✗ dsigmin=1.
67 endif
68 1 call getin('dsigmin', dsigmin)
69 1 WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
70
71
72 ✗ select case (vert_sampling)
73 case ("param")
74 ! On lit les options dans sigma.def:
75 ✗ OPEN(99, file='sigma.def', status='old', form='formatted')
76 ✗ READ(99, *) scaleheight ! hauteur d'echelle 8.
77 ✗ READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
78 ✗ READ(99, *) beta ! facteur d'acroissement en haut 1.3
79 ✗ READ(99, *) k0 ! nombre de couches dans la transition surf
80 ✗ READ(99, *) k1 ! nombre de couches dans la transition haute
81 ✗ CLOSE(99)
82 ✗ alpha=deltaz/(llm*scaleheight)
83 ✗ write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
84 ✗ scaleheight, alpha, k0, k1, beta
85
86 ✗ alpha=deltaz/tanh(1./k0)*2.
87 zkm1=0.
88 ✗ sig(1)=1.
89 ✗ do l=1, llm
90 sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) &
91 *exp(-alpha/scaleheight*tanh((llm-k1)/k0) &
92 ✗ *beta**(l-(llm-k1))/log(beta))
93 ✗ zk=-scaleheight*log(sig(l+1))
94
95 ✗ dzk1=alpha*tanh(l/k0)
96 ✗ dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
97 ✗ write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
98 ✗ zkm1=zk
99 enddo
100
101 ✗ sig(llm+1)=0.
102
103 ✗ bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
104 ✗ bp(llmp1) = 0.
105
106 ✗ ap = pa * (sig - bp)
107 case("sigma")
108 ✗ DO l = 1, llm
109 ✗ x = 2*asin(1.) * (l - 0.5) / (llm + 1)
110 ✗ dsig(l) = dsigmin + 7.0 * SIN(x)**2
111 ENDDO
112 ✗ dsig = dsig / sum(dsig)
113 ✗ sig(llm+1) = 0.
114 ✗ DO l = llm, 1, -1
115 ✗ sig(l) = sig(l+1) + dsig(l)
116 ENDDO
117
118 ✗ bp(1)=1.
119 ✗ bp(2: llm) = sig(2:llm)
120 ✗ bp(llmp1) = 0.
121 ✗ ap(:)=0.
122 case("tropo")
123 ✗ DO l = 1, llm
124 ✗ x = 2*asin(1.) * (l - 0.5) / (llm + 1)
125 ✗ dsig(l) = dsigmin + 7.0 * SIN(x)**2
126 ENDDO
127 ✗ dsig = dsig / sum(dsig)
128 ✗ sig(llm+1) = 0.
129 ✗ DO l = llm, 1, -1
130 ✗ sig(l) = sig(l+1) + dsig(l)
131 ENDDO
132
133 ✗ bp(1)=1.
134 ✗ bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
135 ✗ bp(llmp1) = 0.
136
137 ✗ ap(1)=0.
138 ✗ ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
139 case("strato")
140
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1 times.
 40 DO l = 1, llm
141 39 x = 2*asin(1.) * (l - 0.5) / (llm + 1)
142 dsig(l) =(dsigmin + 7. * SIN(x)**2) &
143 40 *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
144 ENDDO
145
4/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 39 times.
✓ Branch 2 taken 39 times.
✓ Branch 3 taken 1 times.
 80 dsig = dsig / sum(dsig)
146 1 sig(llm+1) = 0.
147
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1 times.
 40 DO l = llm, 1, -1
148 40 sig(l) = sig(l+1) + dsig(l)
149 ENDDO
150
151 1 bp(1)=1.
152
2/2
✓ Branch 0 taken 38 times.
✓ Branch 1 taken 1 times.
 39 bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
153 1 bp(llmp1) = 0.
154
155 1 ap(1)=0.
156
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1 times.
 40 ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
157 case("strato_correct")
158 !==================================================================
159 ! Fredho 2014/05/18, Saint-Louis du Senegal
160 ! Cette version de la discretisation strato est corrige au niveau
161 ! du passage des sig aux ap, bp
162 ! la version precedente donne un coude dans l'epaisseur des couches
163 ! vers la tropopause
164 !==================================================================
165
166
167 ✗ DO l = 1, llm
168 ✗ x = 2*asin(1.) * (l - 0.5) / (llm + 1)
169 dsig(l) =(dsigmin + 7. * SIN(x)**2) &
170 ✗ *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
171 ENDDO
172 ✗ dsig = dsig / sum(dsig)
173 ✗ sig0(llm+1) = 0.
174 ✗ DO l = llm, 1, -1
175 ✗ sig0(l) = sig0(l+1) + dsig(l)
176 ENDDO
177 ✗ sig=racinesig(sig0)
178
179 ✗ bp(1)=1.
180 ✗ bp(2:llm)=EXP(1.-1./sig(2: llm)**2)
181 ✗ bp(llmp1)=0.
182
183 ✗ ap(1)=0.
184 ✗ ap(2:llm)=pa*(sig(2:llm)-bp(2:llm))
185 ✗ ap(llm+1)=0.
186
187 CASE("strato_custom0")
188 !=======================================================
189 ! Version Transitoire
190 ! custumize strato distribution with specific alpha & beta values and function
191 ! depending on llm (experimental and temporary)!
192 SELECT CASE (llm)
193 CASE(55)
194 alpha=0.45
195 beta=4.0
196 CASE(63)
197 alpha=0.45
198 beta=5.0
199 CASE(71)
200 alpha=3.05
201 beta=65.
202 CASE(79)
203 alpha=3.20
204 ! alpha=2.05 ! FLOTT 79 (PLANTE)
205 beta=70.
206 END SELECT
207 ! Or used values provided by user in def file:
208 ✗ CALL getin("strato_alpha",alpha)
209 ✗ CALL getin("strato_beta",beta)
210
211 ! Build geometrical distribution
212 ✗ scaleheight=7.
213 ✗ zz(1)=0.
214 IF (llm==55.OR.llm==63) THEN
215 DO l=1,llm
216 z=zz(l)/scaleheight
217 zz(l+1)=zz(l)+0.03+z*1.5*(1.-TANH(z-0.5))+alpha*(1.+TANH(z-1.5)) &
218 +5.0*EXP((l-llm)/beta)
219 ENDDO
220 ELSEIF (llm==71) THEN !.OR.llm==79) THEN ! FLOTT 79 (PLANTE)
221 DO l=1,llm
222 z=zz(l)
223 zz(l+1)=zz(l)+0.02+0.88*TANH(z/2.5)+alpha*(1.+TANH((z-beta)/15.))
224 ENDDO
225 ELSEIF (llm==79) THEN
226 DO l=1,llm
227 z=zz(l)
228 zz(l+1)=zz(l)+0.02+0.80*TANH(z/3.8)+alpha*(1+TANH((z-beta)/17.)) &
229 +0.03*TANH(z/.25)
230 ENDDO
231 ENDIF ! of IF (llm==55.OR.llm==63) ...
232
233
234 ! Build sigma distribution
235 ✗ sig0=EXP(-zz(:)/scaleheight)
236 ✗ sig0(llm+1)=0.
237 ! sig=ridders(sig0)
238 ✗ sig=racinesig(sig0)
239
240 ! Compute ap() and bp()
241 ✗ bp(1)=1.
242 ✗ bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
243 ✗ bp(llm+1)=0.
244 ✗ ap=pa*(sig-bp)
245
246 CASE("strato_custom")
247 !===================================================================
248 ! David Cugnet, François Lott, Lionel Guez, Ehouoarn Millour, Fredho
249 ! 2014/05
250 ! custumize strato distribution
251 ! Al the parameter are given in km assuming a given scalehigh
252 ✗ vert_scale=7. ! scale hight
253 ✗ vert_dzmin=0.02 ! width of first layer
254 ✗ vert_dzlow=1. ! dz in the low atmosphere
255 ✗ vert_z0low=8. ! height at which resolution recches dzlow
256 ✗ vert_dzmid=3. ! dz in the mid atmsophere
257 ✗ vert_z0mid=70. ! height at which resolution recches dzmid
258 ✗ vert_h_mid=20. ! width of the transition
259 ✗ vert_dzhig=11. ! dz in the high atmsophere
260 ✗ vert_z0hig=80. ! height at which resolution recches dz
261 ✗ vert_h_hig=20. ! width of the transition
262 !===================================================================
263
264 ✗ call getin('vert_scale',vert_scale)
265 ✗ call getin('vert_dzmin',vert_dzmin)
266 ✗ call getin('vert_dzlow',vert_dzlow)
267 ✗ call getin('vert_z0low',vert_z0low)
268 ✗ CALL getin('vert_dzmid',vert_dzmid)
269 ✗ CALL getin('vert_z0mid',vert_z0mid)
270 ✗ call getin('vert_h_mid',vert_h_mid)
271 ✗ call getin('vert_dzhig',vert_dzhig)
272 ✗ call getin('vert_z0hig',vert_z0hig)
273 ✗ call getin('vert_h_hig',vert_h_hig)
274
275 ✗ scaleheight=vert_scale ! for consistency with further computations
276 ! Build geometrical distribution
277 ✗ zz(1)=0.
278 ✗ DO l=1,llm
279 ✗ z=zz(l)
280 zz(l+1)=zz(l)+vert_dzmin+vert_dzlow*TANH(z/vert_z0low)+ &
281 & (vert_dzmid-vert_dzlow)* &
282 & (TANH((z-vert_z0mid)/vert_h_mid)-TANH((-vert_z0mid)/vert_h_mid)) &
283 & +(vert_dzhig-vert_dzmid-vert_dzlow)* &
284 ✗ & (TANH((z-vert_z0hig)/vert_h_hig)-TANH((-vert_z0hig)/vert_h_hig))
285 ENDDO
286
287
288 !===================================================================
289 ! Comment added Fredho 2014/05/18, Saint-Louis, Senegal
290 ! From approximate z to ap, bp, so that p=ap+bp*p0 and p/p0=exp(-z/H)
291 ! sig0 is p/p0
292 ! sig is an intermediate distribution introduce to estimate bp
293 ! 1. sig0=exp(-z/H)
294 ! 2. inversion of sig0=(1-pa/p0)*sig+(1-pa/p0)*exp(1-1/sig**2)
295 ! 3. bp=exp(1-1/sig**2)
296 ! 4. ap deduced from the combination of 2 and 3 so that sig0=ap/p0+bp
297 !===================================================================
298
299 ✗ sig0=EXP(-zz(:)/vert_scale)
300 ✗ sig0(llm+1)=0.
301 ✗ sig=racinesig(sig0)
302 ✗ bp(1)=1.
303 ✗ bp(2:llm)=EXP(1.-1./sig(2:llm)**2)
304 ✗ bp(llm+1)=0.
305 ✗ ap=pa*(sig-bp)
306
307 case("read")
308 ! Read "ap" and "bp". First line is skipped (title line). "ap"
309 ! should be in Pa. First couple of values should correspond to
310 ! the surface, that is : "bp" should be in descending order.
311 ✗ call new_unit(unit)
312 open(unit, file="hybrid.txt", status="old", action="read", &
313 ✗ position="rewind")
314 ✗ read(unit, fmt=*) ! skip title line
315 ✗ do l = 1, llm + 1
316 ✗ read(unit, fmt=*) ap(l), bp(l)
317 end do
318 ✗ close(unit)
319 call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
320 ✗ bp(llm + 1) == 0., "disvert: bad ap or bp values")
321 case default
322
1/9
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1 call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
323 END select
324
325
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1 times.
 40 DO l=1, llm
326 40 nivsigs(l) = REAL(l)
327 ENDDO
328
329
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 1 times.
 41 DO l=1, llmp1
330 41 nivsig(l)= REAL(l)
331 ENDDO
332
333 1 write(lunout, *) trim(modname),': BP '
334 1 write(lunout, *) bp
335 1 write(lunout, *) trim(modname),': AP '
336 1 write(lunout, *) ap
337
338 1 write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
339 1 write(lunout, *)'couches calcules pour une pression de surface =', preff
340 1 write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
341 1 write(lunout, *) scaleheight,' km'
342
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1 times.
 40 DO l = 1, llm
343 39 dpres(l) = bp(l) - bp(l+1)
344 39 aps(l) = 0.5 *( ap(l) +ap(l+1))
345 39 bps(l) = 0.5 *( bp(l) +bp(l+1))
346 39 presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
347 39 pseudoalt(l) = log(preff/presnivs(l))*scaleheight
348 39 write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
349 pseudoalt(l) &
350 39 , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
351 79 max(ap(l+1)+bp(l+1)*preff, 1.e-10))
352 ENDDO
353
354 1 write(lunout, *) trim(modname),': PRESNIVS '
355 1 write(lunout, *) presnivs
356
357 CONTAINS
358
359 !-------------------------------------------------------------------------------
360 !
361 FUNCTION ridders(sig) RESULT(sg)
362 !
363 !-------------------------------------------------------------------------------
364 IMPLICIT NONE
365 !-------------------------------------------------------------------------------
366 ! Purpose: Search for s solving (Pa/Preff)*s+(1-Pa/Preff)*EXP(1-1./s**2)=sg
367 ! Notes: Uses Ridders' method, quite robust. Initial bracketing: 0<=sg<=1.
368 ! Reference: Ridders, C. F. J. "A New Algorithm for Computing a Single Root of a
369 ! Real Continuous Function" IEEE Trans. Circuits Systems 26, 979-980, 1979
370 !-------------------------------------------------------------------------------
371 ! Arguments:
372 REAL, INTENT(IN) :: sig(:)
373 REAL :: sg(SIZE(sig))
374 !-------------------------------------------------------------------------------
375 ! Local variables:
376 INTEGER :: it, ns, maxit
377 REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
378 !-------------------------------------------------------------------------------
379 ns=SIZE(sig); maxit=9999
380 c1=Pa/Preff; c2=1.-c1
381 DO l=1,ns
382 xx=HUGE(1.)
383 x1=0.0; f1=distrib(x1,c1,c2,sig(l))
384 x2=1.0; f2=distrib(x2,c1,c2,sig(l))
385 DO it=1,maxit
386 x3=0.5*(x1+x2); f3=distrib(x3,c1,c2,sig(l))
387 s=SQRT(f3**2-f1*f2); IF(s==0.) EXIT
388 x4=x3+(x3-x1)*(SIGN(1.,f1-f2)*f3/s); IF(ABS(10.*LOG(x4-xx))<=1E-5) EXIT
389 xx=x4; f4=distrib(x4,c1,c2,sig(l)); IF(f4==0.) EXIT
390 IF(SIGN(f3,f4)/=f3) THEN; x1=x3; f1=f3; x2=xx; f2=f4
391 ELSE IF(SIGN(f1,f4)/=f1) THEN; x2=xx; f2=f4
392 ELSE IF(SIGN(f2,f4)/=f2) THEN; x1=xx; f1=f4
393 ELSE; CALL abort_gcm("ridders",'Algorithm failed (which is odd...')
394 END IF
395 IF(ABS(10.*LOG(ABS(x2-x1)))<=1E-5) EXIT !--- ERROR ON SIG <= 0.01m
396 END DO
397 IF(it==maxit+1) WRITE(lunout,'(a,i3)')'WARNING in ridder: failed to converg&
398 &e for level ',l
399 sg(l)=xx
400 END DO
401 sg(1)=1.; sg(ns)=0.
402
403 END FUNCTION ridders
404
405 ✗ FUNCTION racinesig(sig) RESULT(sg)
406 !
407 !-------------------------------------------------------------------------------
408 IMPLICIT NONE
409 !-------------------------------------------------------------------------------
410 ! Fredho 2014/05/18
411 ! Purpose: Search for s solving (Pa/Preff)*sg+(1-Pa/Preff)*EXP(1-1./sg**2)=s
412 ! Notes: Uses Newton Raphson search
413 !-------------------------------------------------------------------------------
414 ! Arguments:
415 REAL, INTENT(IN) :: sig(:)
416 REAL :: sg(SIZE(sig))
417 !-------------------------------------------------------------------------------
418 ! Local variables:
419 INTEGER :: it, ns, maxit
420 REAL :: c1, c2, x1, x2, x3, x4, f1, f2, f3, f4, s, xx, distrib
421 !-------------------------------------------------------------------------------
422 ns=SIZE(sig); maxit=100
423 ✗ c1=Pa/Preff; c2=1.-c1
424 ✗ DO l=2,ns-1
425 ✗ sg(l)=sig(l)
426 ✗ DO it=1,maxit
427 ✗ f1=exp(1-1./sg(l)**2)*(1.-c1)
428 ✗ sg(l)=sg(l)-(c1*sg(l)+f1-sig(l))/(c1+2*f1*sg(l)**(-3))
429 ENDDO
430 ! print*,'SSSSIG ',sig(l),sg(l),c1*sg(l)+exp(1-1./sg(l)**2)*(1.-c1)
431 ENDDO
432 ✗ sg(1)=1.; sg(ns)=0.
433
434 ✗ END FUNCTION racinesig
435
436
437
438
439 END SUBROUTINE disvert
440
441 !-------------------------------------------------------------------------------
442
443 ✗ FUNCTION distrib(x,c1,c2,x0) RESULT(res)
444 !
445 !-------------------------------------------------------------------------------
446 ! Arguments:
447 REAL, INTENT(IN) :: x, c1, c2, x0
448 REAL :: res
449 !-------------------------------------------------------------------------------
450 ✗ res=c1*x+c2*EXP(1-1/(x**2))-x0
451
452 ✗ END FUNCTION distrib
453
454
455