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

Annotation of /trunk/Sources/dyn3d/test_disvert.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 134 - (hide annotations)
Wed Apr 29 15:47:56 2015 UTC (9 years ago) by guez
File size: 3474 byte(s)
Sources inside, compilation outside.
1 guez 99 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 guez 107 use disvert_m, only: ap, bp, preff, presnivs, s
20 guez 99 use dimens_m, only: llm
21     use exner_hyb_m, only: exner_hyb
22     use jumble, only: new_unit
23    
24     integer i, unit, l
25    
26     integer, parameter:: ngrid = 8
27    
28     real p(ngrid, 1, llm + 1) ! pressure at half-level, in Pa
29     real z(llm) ! pressure-altitude at half-level (km)
30     real pks(ngrid, 1) ! exner function at the surface, in J K-1 kg-1
31     real pk(ngrid, 1, llm) ! exner function at full level, in J K-1 kg-1
32     real ps(ngrid, 1) ! surface pressure, in Pa
33     real p_lay(ngrid, llm) ! pressure at full level, in Pa
34     real, parameter:: delta_ps = 6e4 / (ngrid - 2) ! in Pa
35    
36     !---------------------
37    
38     print *, "Call sequence information: test_disvert"
39    
40     ps(:, 1) = (/(5e4 + delta_ps * i, i = 0, ngrid - 2), preff/)
41     forall (l = 1: llm + 1) p(:, 1, l) = ap(l) + bp(l) * ps(:, 1)
42     call exner_hyb(ps, p, pks, pk)
43     p_lay = preff * (pk(:, 1, :) / cpp)**(1. / kappa)
44    
45     ! Write distribution for the reference surface pressure (index ngrid):
46    
47     z = 7. * log(preff / p(ngrid, 1, :llm))
48     call new_unit(unit)
49    
50     open(unit, file="half_level.csv", status="replace", action="write")
51     ! Title line:
52     write(unit, fmt=*) '"ap (hPa)" "bp" "s" "pressure (hPa)" "z (km)"'
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 line:
61     write(unit, fmt=*) &
62     '"pressure (hPa)" "z (km)" "presnivs (hPa)" "delta z (km)"'
63     do l = 1, llm - 1
64     write(unit, fmt=*) p_lay(ngrid, l) / 100., &
65     7. * log(preff / p_lay(ngrid, l)), presnivs(l) / 100., z(l+1) - z(l)
66     end do
67     write(unit, fmt=*) p_lay(ngrid, llm) / 100., &
68     7. * log(preff / p_lay(ngrid, llm)), presnivs(llm) / 100.
69     close(unit)
70     print *, 'The file "full_level.csv" has been created.'
71    
72 guez 116 ! Are pressure values in the right order?
73     if (any(p(:, 1, :llm) <= p_lay .or. p_lay <= p(:, 1, 2:))) then
74     ! List details and stop:
75     do l = 1, llm
76     do i = 1, ngrid
77     if (p(i, 1, l) <= p_lay(i, l)) then
78     print 1000, "ps = ", ps(i, 1) / 100., "hPa, p(level ", l, &
79     ") = ", p(i, 1, l) / 100., " hPa <= p(layer ", l, ") = ", &
80     p_lay(i, l) / 100., " hPa"
81     end if
82     if (p_lay(i, l) <= p(i, 1, l+1)) then
83     print 1000, "ps = ", ps(i, 1) / 100., &
84     "hPa, p(layer ", l, ") = ", p_lay(i, l) / 100., &
85     " hPa <= p(level ", l + 1, ") = ", &
86     p(i, 1, l + 1) / 100., " hPa"
87     end if
88     end do
89     end do
90     call abort_gcm("test_disvert", "bad order of pressure values", 1)
91     end if
92    
93 guez 99 1000 format (3(a, g10.4, a, i0))
94    
95     end subroutine test_disvert
96    
97     end module test_disvert_m

  ViewVC Help
Powered by ViewVC 1.1.21