trunk/dyn3d/disvert.f

module disvert_m

  use dimens_m, only: llm

  implicit none

  private llm

  real ap(llm+1), pa ! in Pa
  real bp(llm+1)
  real presnivs(llm) ! pressions approximatives des milieux de couches, en Pa
  real, parameter:: preff = 101325. ! in Pa
  real nivsigs(llm), nivsig(llm+1)

  save

contains

  SUBROUTINE disvert

    ! From dyn3d/disvert.F, v 1.1.1.1 2004/05/19 12:53:05
    ! Author: P. Le Van

    ! This procedure sets the vertical grid. It defines the host
    ! variables "ap", "bp", "presnivs", "nivsigs" and "nivsig". "pa"
    ! should be defined before this procedure is called.

    use jumble, only: new_unit
    use nr_util, only: pi, assert
    use unit_nml_m, only: unit_nml

    ! Local:

    REAL s(llm+1)
    ! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
    ! the interface between layers "l" and "l-1"

    real ds(llm)
    ! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"

    INTEGER l, unit
    REAL alpha, x(llm)

    character(len=7):: vert_sampling = "tropo"
    ! other allowed values are "param", "strato" and "read"

    real:: h = 7. ! scale height, in km
    ! used only if vert_sampling == "param"

    ! These variables are used only in the case vert_sampling == "param":
    real:: deltaz = 0.04 ! épaisseur de la première couche
    real:: beta = 1.3 ! facteur d'accroissement en haut
    real:: k0 = 20. ! nombre de couches dans la transition surface
    real:: k1 = 1.2 ! nombre de couches dans la transition haute

    namelist /disvert_nml/h, deltaz, beta, k0, k1, vert_sampling

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

    print *, "Call sequence information: disvert"

    forall (l = 1: llm) nivsigs(l) = REAL(l)
    forall (l = 1: llm + 1) nivsig(l) = REAL(l)

    print *, "Enter namelist 'disvert_nml'."
    read(unit=*, nml=disvert_nml)
    write(unit_nml, nml=disvert_nml)

    select case (vert_sampling)
    case ("param")
       s(1) = 1.
       s(llm+1) = 0.
       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))
       call compute_ab
    case ("tropo")
       s(1) = 1.
       s(llm+1) = 0.
       forall (l = 1: llm) ds(l) &
            = 1. + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2
       ds = ds / sum(ds)

       DO l = llm, 2, -1
          s(l) = s(l+1) + ds(l)
       ENDDO

       call compute_ab
    case ("strato")
       ! Recommended by F. Lott for a domain including the stratosphere
       s(1) = 1.
       s(llm+1) = 0.
       forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)

       ds = (0.3 + 7. * SIN(x)**2) * (1. - tanh(2 * x / pi - 1.))**2 / 4.
       ds = ds / sum(ds)

       DO l = llm, 2, -1
          s(l) = s(l+1) + ds(l)
       ENDDO

       call compute_ab
    case("read")
       ! Read "ap" and "bp". First line is skipped (title line). "ap"
       ! should be in Pa. First couple of values should correspond to
       ! the surface, that is : "bp" should be in descending order.
       call new_unit(unit)
       open(unit, file="hybrid.csv", status="old", action="read", &
            position="rewind")
       read(unit, fmt=*) ! skip title line
       do l = 1, llm + 1
          read(unit, fmt=*) ap(l), bp(l)
       end do
       close(unit)
       ! Quick check:
       call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
            bp(llm + 1) == 0., "disvert: bad ap or bp values")
    case default
       print *, 'Wrong value for "vert_sampling"'
       stop 1
    END select

    forall (l = 1: llm) presnivs(l) = 0.5 &
         * (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)

  contains

    subroutine compute_ab

      ! Calcul de "ap" et "bp" :
      bp(:llm) = EXP(1. - 1. / s(:llm)**2)
      bp(llm + 1) = 0.
      ap = pa * (s - bp)

    end subroutine compute_ab

  END SUBROUTINE disvert

end module disvert_m
1	guez	66	module disvert_m
2	guez	3
3			use dimens_m, only: llm
4
5			implicit none
6
7			private llm
8
9	guez	20	real ap(llm+1), pa ! in Pa
10	guez	32	real bp(llm+1)
11	guez	22	real presnivs(llm) ! pressions approximatives des milieux de couches, en Pa
12	guez	3	real, parameter:: preff = 101325. ! in Pa
13			real nivsigs(llm), nivsig(llm+1)
14
15			save
16
17			contains
18
19			SUBROUTINE disvert
20
21			! From dyn3d/disvert.F, v 1.1.1.1 2004/05/19 12:53:05
22	guez	66	! Author: P. Le Van
23	guez	3
24	guez	66	! This procedure sets the vertical grid. It defines the host
25			! variables "ap", "bp", "presnivs", "nivsigs" and "nivsig". "pa"
26			! should be defined before this procedure is called.
27	guez	3
28	guez	66	use jumble, only: new_unit
29	guez	78	use nr_util, only: pi, assert
30	guez	57	use unit_nml_m, only: unit_nml
31	guez	3
32	guez	66	! Local:
33
34	guez	3	REAL s(llm+1)
35	guez	22	! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
36			! the interface between layers "l" and "l-1"
37	guez	3
38			real ds(llm)
39			! "ds(l)" : épaisseur de la couche "l" dans la coordonnée "s"
40
41	guez	42	INTEGER l, unit
42	guez	53	REAL alpha, x(llm)
43	guez	3
44	guez	78	character(len=7):: vert_sampling = "tropo"
45			! other allowed values are "param", "strato" and "read"
46	guez	3
47			real:: h = 7. ! scale height, in km
48	guez	78	! used only if vert_sampling == "param"
49	guez	3
50	guez	78	! These variables are used only in the case vert_sampling == "param":
51	guez	3	real:: deltaz = 0.04 ! épaisseur de la première couche
52			real:: beta = 1.3 ! facteur d'accroissement en haut
53			real:: k0 = 20. ! nombre de couches dans la transition surface
54			real:: k1 = 1.2 ! nombre de couches dans la transition haute
55
56	guez	78	namelist /disvert_nml/h, deltaz, beta, k0, k1, vert_sampling
57	guez	3
58			!-----------------------------------------------------------------------
59
60			print *, "Call sequence information: disvert"
61
62			forall (l = 1: llm) nivsigs(l) = REAL(l)
63			forall (l = 1: llm + 1) nivsig(l) = REAL(l)
64
65			print *, "Enter namelist 'disvert_nml'."
66			read(unit=*, nml=disvert_nml)
67	guez	57	write(unit_nml, nml=disvert_nml)
68	guez	3
69	guez	78	select case (vert_sampling)
70	guez	3	case ("param")
71	guez	42	s(1) = 1.
72			s(llm+1) = 0.
73	guez	3	alpha = deltaz / tanh(1./k0) * 2.
74			forall (l = 2: llm) s(l) &
75			= cosh((l - 1) / k0) *(- alpha k0 / h) &
76			* exp(- alpha / h * tanh((llm - k1) / k0) &
77			* beta **(l - 1 - (llm - k1)) / log(beta))
78	guez	78	call compute_ab
79	guez	65	case ("tropo")
80	guez	42	s(1) = 1.
81			s(llm+1) = 0.
82	guez	3	forall (l = 1: llm) ds(l) &
83	guez	66	= 1. + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2
84	guez	3	ds = ds / sum(ds)
85
86			DO l = llm, 2, -1
87			s(l) = s(l+1) + ds(l)
88			ENDDO
89
90	guez	78	call compute_ab
91	guez	65	case ("strato")
92	guez	3	! Recommended by F. Lott for a domain including the stratosphere
93	guez	42	s(1) = 1.
94			s(llm+1) = 0.
95	guez	3	forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
96
97	guez	66	ds = (0.3 + 7. * SIN(x)*2) (1. - tanh(2 * x / pi - 1.))**2 / 4.
98	guez	3	ds = ds / sum(ds)
99
100			DO l = llm, 2, -1
101			s(l) = s(l+1) + ds(l)
102			ENDDO
103	guez	78
104			call compute_ab
105	guez	42	case("read")
106	guez	78	! Read "ap" and "bp". First line is skipped (title line). "ap"
107			! should be in Pa. First couple of values should correspond to
108			! the surface, that is : "bp" should be in descending order.
109	guez	42	call new_unit(unit)
110			open(unit, file="hybrid.csv", status="old", action="read", &
111			position="rewind")
112			read(unit, fmt=*) ! skip title line
113			do l = 1, llm + 1
114	guez	78	read(unit, fmt=*) ap(l), bp(l)
115	guez	42	end do
116			close(unit)
117			! Quick check:
118	guez	78	call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
119			bp(llm + 1) == 0., "disvert: bad ap or bp values")
120	guez	3	case default
121	guez	78	print *, 'Wrong value for "vert_sampling"'
122	guez	3	stop 1
123			END select
124
125	guez	32	forall (l = 1: llm) presnivs(l) = 0.5 &
126			* (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
127	guez	3
128	guez	78	contains
129
130			subroutine compute_ab
131
132			! Calcul de "ap" et "bp" :
133			bp(:llm) = EXP(1. - 1. / s(:llm)**2)
134			bp(llm + 1) = 0.
135			ap = pa * (s - bp)
136
137			end subroutine compute_ab
138
139	guez	3	END SUBROUTINE disvert
140
141	guez	66	end module disvert_m