trunk/dyn3d/test_disvert.f

module test_disvert_m

  implicit none

contains

  subroutine test_disvert

    ! Author: Lionel GUEZ

    ! This procedure tests the order of pressure values at half-levels
    ! and full levels, and writes the vertical coordinates to files
    ! "half_level.csv" and "full_level.csv" for a reference surface
    ! pressure. We arbitrarily choose to test ngrid values of the
    ! surface pressure, which sample possible values on Earth.

    use abort_gcm_m, only: abort_gcm
    use comconst, only: kappa, cpp
    use disvert_m, only: ap, bp, preff, presnivs, s
    use dimens_m, only: llm
    use exner_hyb_m, only: exner_hyb
    use jumble, only: new_unit

    integer i, unit, l

    integer, parameter:: ngrid = 8

    real p(ngrid, 1, llm + 1) ! pressure at half-level, in Pa
    real z(llm) ! pressure-altitude at half-level (km)
    real pks(ngrid, 1) ! exner function at the surface, in J K-1 kg-1 
    real pk(ngrid, 1, llm) ! exner function at full level, in J K-1 kg-1 
    real ps(ngrid, 1) ! surface pressure, in Pa
    real p_lay(ngrid, llm) ! pressure at full level, in Pa
    real, parameter:: delta_ps = 6e4 / (ngrid - 2) ! in Pa

    !---------------------

    print *, "Call sequence information: test_disvert"

    ps(:, 1) = (/(5e4 + delta_ps * i, i = 0, ngrid - 2), preff/)
    forall (l = 1: llm + 1) p(:, 1, l) = ap(l) + bp(l) * ps(:, 1)
    call exner_hyb(ps, p, pks, pk)
    p_lay = preff * (pk(:, 1, :) / cpp)**(1. / kappa)

    ! Write distribution for the reference surface pressure (index ngrid):

    z = 7. * log(preff / p(ngrid, 1, :llm))
    call new_unit(unit)

    open(unit, file="half_level.csv", status="replace", action="write")
    ! Title line:
    write(unit, fmt=*) '"ap (hPa)" "bp" "s" "pressure (hPa)" "z (km)"'
    do l = 1, llm
       write(unit, fmt=*) ap(l) / 100., bp(l), s(l), p(ngrid, 1, l) / 100., z(l)
    end do
    close(unit)
    print *, 'The file "half_level.csv" has been created.'

    open(unit, file="full_level.csv", status="replace", action="write")
    ! Title line:
    write(unit, fmt=*) &
         '"pressure (hPa)" "z (km)" "presnivs (hPa)" "delta z (km)"'
    do l = 1, llm - 1
       write(unit, fmt=*) p_lay(ngrid, l) / 100., &
            7. * log(preff / p_lay(ngrid, l)), presnivs(l) / 100., z(l+1) - z(l)
    end do
    write(unit, fmt=*) p_lay(ngrid, llm) / 100., &
         7. * log(preff / p_lay(ngrid, llm)), presnivs(llm) / 100.
    close(unit)
    print *, 'The file "full_level.csv" has been created.'

    ! Are pressure values in the right order?
    if (any(p(:, 1, :llm) <= p_lay .or. p_lay <= p(:, 1, 2:))) then
       ! List details and stop:
       do l = 1, llm
          do i = 1, ngrid
             if (p(i, 1, l) <= p_lay(i, l)) then
                print 1000, "ps = ", ps(i, 1) / 100., "hPa, p(level ",  l, &
                     ") = ", p(i, 1, l) / 100., " hPa <= p(layer ", l, ") = ", &
                     p_lay(i, l) / 100., " hPa"
             end if
             if (p_lay(i, l) <= p(i, 1, l+1)) then
                print 1000, "ps = ", ps(i, 1) / 100., &
                     "hPa, p(layer ", l, ") = ", p_lay(i, l) / 100., &
                     " hPa <= p(level ", l + 1, ") = ", &
                     p(i, 1, l + 1) / 100., " hPa"
             end if
          end do
       end do
       call abort_gcm("test_disvert", "bad order of pressure values")
    end if

1000 format (3(a, g11.4, a, i0))

  end subroutine test_disvert

end module test_disvert_m
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	guez	171	call abort_gcm("test_disvert", "bad order of pressure values")
91	guez	116	end if
92
93	guez	175	1000 format (3(a, g11.4, a, i0))
94	guez	99
95			end subroutine test_disvert
96
97			end module test_disvert_m