/[lmdze]/trunk/dyn3d/test_disvert.f
ViewVC logotype

Contents of /trunk/dyn3d/test_disvert.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 266 - (show annotations)
Thu Apr 19 17:54:55 2018 UTC (6 years ago) by guez
File size: 3753 byte(s)
Define macros of the preprocessor CPP_IIM, CPP_JJM, CPP_LLM so we can
control the resolution from the compilation command, and automate
compilation for several resolutions.

In module yoethf_m, transform variables into named constants. So we do
not need procedure yoethf any longer.

Bug fix in program test_inter_barxy, missing calls to fyhyp and fxhyp,
and definition of rlatu.

Remove variable iecri of module conf_gcm_m. The files dyn_hist*.nc are
written every time step. We are simplifying the output system, pending
replacement by a whole new system.

Modify possible value of vert_sampling from "param" to
"strato_custom", following LMDZ. Default values of corresponding
namelist variables are now the values used for LMDZ CMIP6.

1 module test_disvert_m
2
3 implicit none
4
5 contains
6
7 subroutine test_disvert
8
9 ! Author: Lionel GUEZ
10
11 ! This procedure tests the order of pressure values at half-levels
12 ! and full levels, and writes the vertical coordinates to files
13 ! "half_level.csv" and "full_level.csv" for a reference surface
14 ! pressure. We arbitrarily choose to test ngrid values of the
15 ! surface pressure, which sample possible values on Earth.
16
17 use abort_gcm_m, only: abort_gcm
18 use comconst, only: kappa, cpp
19 use disvert_m, only: ap, bp, preff, presnivs, s
20 use dimensions, only: llm
21 use exner_hyb_m, only: exner_hyb
22 use jumble, only: new_unit
23
24 integer i, unit, l
25 integer, parameter:: ngrid = 8
26 real p(ngrid, 1, llm + 1) ! pressure at half-level, in Pa
27 real z(llm) ! pressure-altitude at half-level (km)
28 real pks(ngrid, 1) ! exner function at the surface, in J K-1 kg-1
29 real pk(ngrid, 1, llm) ! exner function at full level, in J K-1 kg-1
30 real ps(ngrid, 1) ! surface pressure, in Pa
31 real p_lay(ngrid, llm) ! pressure at full level, in Pa
32 real, parameter:: delta_ps = 6e4 / (ngrid - 2) ! in Pa
33
34 !---------------------
35
36 print *, "Call sequence information: test_disvert"
37
38 ps(:, 1) = [(5e4 + delta_ps * i, i = 0, ngrid - 2), preff]
39 forall (l = 1: llm + 1) p(:, 1, l) = ap(l) + bp(l) * ps(:, 1)
40 call exner_hyb(ps, p, pks, pk)
41 p_lay = preff * (pk(:, 1, :) / cpp)**(1. / kappa)
42
43 ! Write distribution for the reference surface pressure (which is
44 ! surface pressure at index ngrid):
45
46 z = 7. * log(preff / p(ngrid, 1, :llm))
47 call new_unit(unit)
48
49 open(unit, file="half_level.csv", status="replace", action="write")
50 ! Title lines:
51 write(unit, fmt=*) 'hPa "" "" hPa km'
52 write(unit, fmt=*) 'ap bp s pressure z'
53 do l = 1, llm
54 write(unit, fmt=*) ap(l) / 100., bp(l), s(l), p(ngrid, 1, l) / 100., z(l)
55 end do
56 close(unit)
57 print *, 'The file "half_level.csv" has been created.'
58
59 open(unit, file="full_level.csv", status="replace", action="write")
60 ! Title lines:
61 write(unit, fmt=*) '"pressure at full level from Exner function for a ', &
62 'reference surface pressure" altitude "approximate pressure at ', &
63 'full level for a reference surface pressure" "layer thickness"'
64 write(unit, fmt=*) 'hPa km hPa km'
65 write(unit, fmt=*) 'p_lay z presnivs delta_z'
66 do l = 1, llm - 1
67 write(unit, fmt=*) p_lay(ngrid, l) / 100., &
68 7. * log(preff / p_lay(ngrid, l)), presnivs(l) / 100., z(l+1) - z(l)
69 end do
70 write(unit, fmt=*) p_lay(ngrid, llm) / 100., &
71 7. * log(preff / p_lay(ngrid, llm)), presnivs(llm) / 100., "NaN"
72 close(unit)
73 print *, 'The file "full_level.csv" has been created.'
74
75 ! Are pressure values in the right order?
76 if (any(p(:, 1, :llm) <= p_lay .or. p_lay <= p(:, 1, 2:))) then
77 ! List details and stop:
78 do l = 1, llm
79 do i = 1, ngrid
80 if (p(i, 1, l) <= p_lay(i, l)) then
81 print 1000, "ps = ", ps(i, 1) / 100., "hPa, p(level ", l, &
82 ") = ", p(i, 1, l) / 100., " hPa <= p(layer ", l, ") = ", &
83 p_lay(i, l) / 100., " hPa"
84 end if
85 if (p_lay(i, l) <= p(i, 1, l+1)) then
86 print 1000, "ps = ", ps(i, 1) / 100., &
87 "hPa, p(layer ", l, ") = ", p_lay(i, l) / 100., &
88 " hPa <= p(level ", l + 1, ") = ", &
89 p(i, 1, l + 1) / 100., " hPa"
90 end if
91 end do
92 end do
93 call abort_gcm("test_disvert", "bad order of pressure values")
94 end if
95
96 1000 format (3(a, g11.4, a, i0))
97
98 end subroutine test_disvert
99
100 end module test_disvert_m

  ViewVC Help
Powered by ViewVC 1.1.21