trunk/dyn3d/disvert.f

module disvert_m

  use dimensions, only: llm

  implicit none

  private llm, hybrid, funcd, y, ya, compute_ab

  real, save:: ap(llm+1) ! in Pa
  real, save:: bp(llm+1)

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

  real, save:: presnivs(llm)
  ! approximate full level pressure for a reference surface pressure, in Pa

  real, parameter:: preff = 101325. ! in Pa
  real y, ya ! for the hybrid function

contains

  SUBROUTINE disvert

    ! From dyn3d/disvert.F, version 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". "pa" should be defined before
    ! this procedure is called.

    ! Libraries:
    use jumble, only: read_column, new_unit
    use nr_util, only: pi, assert
    
    use dynetat0_chosen_m, only: pa
    use unit_nml_m, only: unit_nml

    ! Local:

    real ds(llm)
    ! "ds(l)" : \'epaisseur de la couche "l" dans la coordonn\'ee "s"

    INTEGER l, unit
    REAL x(llm)
    real:: dsigmin = 1.
    real zz(llm) ! in km

    character(len=20):: vert_sampling = "tropo"
    ! Allowed values: "tropo", "strato_custom", "strato",
    ! "read_hybrid", "read_pressure".

    ! These variables are used only in the case vert_sampling ==
    ! "strato_custom", and all are in km:
    real:: vert_scale = 7. ! scale height
    real:: vert_dzmin = 0.017 ! width of first layer
    real:: vert_dzlow = 1. ! dz in the low atmosphere
    real:: vert_z0low = 8.7 ! height at which resolution reaches dzlow
    real:: vert_dzmid = 2. ! dz in the mid atmosphere
    real:: vert_z0mid = 70. ! height at which resolution reaches dzmid
    real:: vert_h_mid = 20. ! width of the transition
    real:: vert_dzhig = 11. ! dz in the high atmosphere
    real:: vert_z0hig = 75. ! height at which resolution reaches dz
    real:: vert_h_hig = 20. ! width of the transition

    real, allocatable:: p(:) ! (2:llm or llm + 1) pressure (in hPa)

    namelist /disvert_nml/vert_sampling, vert_scale, vert_dzmin, vert_dzlow, &
         vert_z0low, vert_dzmid, vert_z0mid, vert_h_mid, vert_dzhig, &
         vert_z0hig, vert_h_hig, dsigmin

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

    print *, "Call sequence information: disvert"

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

    s(1) = 1.
    s(llm+1) = 0.

    select case (vert_sampling)

    case ("tropo")
       ! with llm = 19 and dsigmin = 1 for CMIP 3

       forall (l = 1: llm) ds(l) &
            = dsigmin + 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")
       ! with llm = 39 and dsigmin = 0.3 for CMIP5

       forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)

       ds = (dsigmin + 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 ("strato_custom")
       ! with llm = 79 for CMIP 6

       zz(1) = 0.
       DO l = 1, llm - 1
          zz(l + 1) = zz(l) + vert_dzmin + vert_dzlow &
               * TANH(zz(l) / vert_z0low) + (vert_dzmid - vert_dzlow) &
               * (TANH((zz(l) - vert_z0mid) / vert_h_mid) &
               - TANH(- vert_z0mid / vert_h_mid)) &
               + (vert_dzhig - vert_dzmid - vert_dzlow) &
               * (TANH((zz(l) - vert_z0hig) / vert_h_hig) &
               - TANH(- vert_z0hig / vert_h_hig))
       ENDDO

       allocate(p(2: llm))
       p = preff * EXP(- zz(2:) / vert_scale)
       ya = pa / preff
       s(2: llm) = hybrid(p)

       call compute_ab

    case("read_hybrid")
       ! 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")
       s(2: llm) = ap(2: llm) / pa + bp(2: llm)

    case("read_pressure")
       ! Read pressure values, in Pa, in descending order, from preff
       ! to 0. First line is skipped (title line).
       call read_column("pressure.txt", p, first=2)
       call assert(size(p) == llm + 1, "disvert: bad number of pressure values")
       ! Quick check:
       call assert(p(1) == preff, p(llm + 1) == 0., &
            "disvert: bad pressure values")
       ya = pa / preff
       s(2: llm) = hybrid(p(2: llm))
       call compute_ab

    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)

  END SUBROUTINE disvert

  !**********************************************************

  subroutine compute_ab

    use dynetat0_chosen_m, only: pa

    ! Calcul de "ap" et "bp".

    where (s >= 1. / sqrt(1. - log(tiny(0.))))
       bp = exp(1. - 1. / s**2)
    elsewhere
       bp = 0.
    end where

    ap = pa * (s - bp)

  end subroutine compute_ab

  !**********************************************************

  function hybrid(p)

    ! This procedure computes the hybrid sigma-pressure coordinate
    ! from pressure values, assuming some reference surface
    ! pressure. The procedure assumes, and does not check, that
    ! pressure values are given in descending order.

    use numer_rec_95, only: rtsafe

    real, intent(in):: p(:) ! pressure (in hPa)
    real hybrid(size(p)) ! hybrid sigma-pressure coordinate

    ! Local:
    integer l

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

    y = p(1) / preff
    hybrid(1) = rtsafe(funcd, x1 = 0., x2 = 1., xacc = 1e-4)

    do l = 2, size(p)
       y = p(l) / preff
       ! Assuming descending order in pressure:
       hybrid(l) = rtsafe(funcd, x1 = 0., x2 = hybrid(l - 1), &
            xacc = hybrid(l - 1) * 1e-4)
    end do

  end function hybrid

  !**********************************************************

  SUBROUTINE funcd(s, fval, fderiv)

    REAL, INTENT(IN):: s
    REAL, INTENT(OUT):: fval, fderiv

    ! Local:
    real b

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

    if (s**3 > 1. / huge(0.)) then
       b = exp(1. - 1. / s**2)
       fval = ya * s + (1. - ya) * b - y
       fderiv = ya + 2 * (1. - ya) * b / s**3
    else
       fval = ya * s - y
       fderiv = ya
    end if

  END SUBROUTINE funcd

end module disvert_m
1	guez	66	module disvert_m
2	guez	3
3	guez	265	use dimensions, only: llm
4	guez	3
5			implicit none
6
7	guez	99	private llm, hybrid, funcd, y, ya, compute_ab
8	guez	3
9	guez	313	real, save:: ap(llm+1) ! in Pa
10	guez	83	real, save:: bp(llm+1)
11
12	guez	99	REAL s(llm+1)
13			! "s(l)" is the atmospheric hybrid sigma-pressure coordinate at
14			! half-level, between layers "l" and "l-1"
15
16	guez	83	real, save:: presnivs(llm)
17	guez	266	! approximate full level pressure for a reference surface pressure, in Pa
18	guez	83
19	guez	3	real, parameter:: preff = 101325. ! in Pa
20	guez	99	real y, ya ! for the hybrid function
21	guez	3
22			contains
23
24			SUBROUTINE disvert
25
26	guez	90	! From dyn3d/disvert.F, version 1.1.1.1, 2004/05/19 12:53:05
27	guez	66	! Author: P. Le Van
28	guez	3
29	guez	66	! This procedure sets the vertical grid. It defines the host
30	guez	83	! variables "ap", "bp", "presnivs". "pa" should be defined before
31			! this procedure is called.
32	guez	3
33	guez	313	! Libraries:
34	guez	99	use jumble, only: read_column, new_unit
35	guez	78	use nr_util, only: pi, assert
36	guez	313
37			use dynetat0_chosen_m, only: pa
38	guez	57	use unit_nml_m, only: unit_nml
39	guez	3
40	guez	66	! Local:
41
42	guez	3	real ds(llm)
43	guez	90	! "ds(l)" : \'epaisseur de la couche "l" dans la coordonn\'ee "s"
44	guez	3
45	guez	42	INTEGER l, unit
46	guez	99	REAL x(llm)
47			real:: dsigmin = 1.
48			real zz(llm) ! in km
49	guez	3
50	guez	99	character(len=20):: vert_sampling = "tropo"
51	guez	266	! Allowed values: "tropo", "strato_custom", "strato",
52			! "read_hybrid", "read_pressure".
53	guez	3
54	guez	99	! These variables are used only in the case vert_sampling ==
55	guez	266	! "strato_custom", and all are in km:
56	guez	99	real:: vert_scale = 7. ! scale height
57	guez	266	real:: vert_dzmin = 0.017 ! width of first layer
58	guez	99	real:: vert_dzlow = 1. ! dz in the low atmosphere
59	guez	266	real:: vert_z0low = 8.7 ! height at which resolution reaches dzlow
60			real:: vert_dzmid = 2. ! dz in the mid atmosphere
61	guez	99	real:: vert_z0mid = 70. ! height at which resolution reaches dzmid
62			real:: vert_h_mid = 20. ! width of the transition
63			real:: vert_dzhig = 11. ! dz in the high atmosphere
64	guez	266	real:: vert_z0hig = 75. ! height at which resolution reaches dz
65	guez	99	real:: vert_h_hig = 20. ! width of the transition
66	guez	3
67	guez	266	real, allocatable:: p(:) ! (2:llm or llm + 1) pressure (in hPa)
68	guez	3
69	guez	99	namelist /disvert_nml/vert_sampling, vert_scale, vert_dzmin, vert_dzlow, &
70			vert_z0low, vert_dzmid, vert_z0mid, vert_h_mid, vert_dzhig, &
71			vert_z0hig, vert_h_hig, dsigmin
72	guez	3
73			!-----------------------------------------------------------------------
74
75			print *, "Call sequence information: disvert"
76
77	guez	99	write(unit=*, nml=disvert_nml)
78	guez	3	print *, "Enter namelist 'disvert_nml'."
79			read(unit=*, nml=disvert_nml)
80	guez	57	write(unit_nml, nml=disvert_nml)
81	guez	3
82	guez	99	s(1) = 1.
83			s(llm+1) = 0.
84
85	guez	78	select case (vert_sampling)
86	guez	99
87	guez	65	case ("tropo")
88	guez	266	! with llm = 19 and dsigmin = 1 for CMIP 3
89	guez	99
90	guez	3	forall (l = 1: llm) ds(l) &
91	guez	99	= dsigmin + 7. * SIN(pi * (REAL(l) - 0.5) / real(llm + 1))**2
92	guez	3	ds = ds / sum(ds)
93
94			DO l = llm, 2, -1
95			s(l) = s(l+1) + ds(l)
96			ENDDO
97
98	guez	78	call compute_ab
99	guez	99
100	guez	65	case ("strato")
101	guez	266	! with llm = 39 and dsigmin = 0.3 for CMIP5
102	guez	99
103	guez	3	forall (l = 1: llm) x(l) = pi * (l - 0.5) / (llm + 1)
104
105	guez	99	ds = (dsigmin + 7. * SIN(x)*2) (1. - tanh(2 * x / pi - 1.))**2 / 4.
106	guez	3	ds = ds / sum(ds)
107
108			DO l = llm, 2, -1
109			s(l) = s(l+1) + ds(l)
110			ENDDO
111	guez	78
112			call compute_ab
113	guez	99
114	guez	266	case ("strato_custom")
115	guez	99	! with llm = 79 for CMIP 6
116
117			zz(1) = 0.
118			DO l = 1, llm - 1
119			zz(l + 1) = zz(l) + vert_dzmin + vert_dzlow &
120			* TANH(zz(l) / vert_z0low) + (vert_dzmid - vert_dzlow) &
121			* (TANH((zz(l) - vert_z0mid) / vert_h_mid) &
122			- TANH(- vert_z0mid / vert_h_mid)) &
123			+ (vert_dzhig - vert_dzmid - vert_dzlow) &
124			* (TANH((zz(l) - vert_z0hig) / vert_h_hig) &
125			- TANH(- vert_z0hig / vert_h_hig))
126			ENDDO
127
128			allocate(p(2: llm))
129			p = preff * EXP(- zz(2:) / vert_scale)
130			ya = pa / preff
131			s(2: llm) = hybrid(p)
132	guez	253
133	guez	99	call compute_ab
134
135			case("read_hybrid")
136	guez	78	! Read "ap" and "bp". First line is skipped (title line). "ap"
137			! should be in Pa. First couple of values should correspond to
138			! the surface, that is : "bp" should be in descending order.
139	guez	42	call new_unit(unit)
140			open(unit, file="hybrid.csv", status="old", action="read", &
141			position="rewind")
142			read(unit, fmt=*) ! skip title line
143			do l = 1, llm + 1
144	guez	78	read(unit, fmt=*) ap(l), bp(l)
145	guez	42	end do
146			close(unit)
147			! Quick check:
148	guez	78	call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
149			bp(llm + 1) == 0., "disvert: bad ap or bp values")
150	guez	99	s(2: llm) = ap(2: llm) / pa + bp(2: llm)
151
152			case("read_pressure")
153			! Read pressure values, in Pa, in descending order, from preff
154			! to 0. First line is skipped (title line).
155			call read_column("pressure.txt", p, first=2)
156			call assert(size(p) == llm + 1, "disvert: bad number of pressure values")
157			! Quick check:
158			call assert(p(1) == preff, p(llm + 1) == 0., &
159			"disvert: bad pressure values")
160			ya = pa / preff
161			s(2: llm) = hybrid(p(2: llm))
162			call compute_ab
163
164	guez	3	case default
165	guez	78	print *, 'Wrong value for "vert_sampling"'
166	guez	3	stop 1
167	guez	99
168	guez	3	END select
169
170	guez	32	forall (l = 1: llm) presnivs(l) = 0.5 &
171			* (ap(l) + bp(l) * preff + ap(l+1) + bp(l+1) * preff)
172	guez	3
173	guez	99	END SUBROUTINE disvert
174	guez	78
175	guez	99	!**********************************************************
176	guez	78
177	guez	99	subroutine compute_ab
178	guez	78
179	guez	313	use dynetat0_chosen_m, only: pa
180
181	guez	99	! Calcul de "ap" et "bp".
182	guez	78
183	guez	99	where (s >= 1. / sqrt(1. - log(tiny(0.))))
184			bp = exp(1. - 1. / s**2)
185			elsewhere
186			bp = 0.
187			end where
188	guez	3
189	guez	99	ap = pa * (s - bp)
190
191			end subroutine compute_ab
192
193			!**********************************************************
194
195			function hybrid(p)
196
197			! This procedure computes the hybrid sigma-pressure coordinate
198			! from pressure values, assuming some reference surface
199			! pressure. The procedure assumes, and does not check, that
200			! pressure values are given in descending order.
201
202			use numer_rec_95, only: rtsafe
203
204			real, intent(in):: p(:) ! pressure (in hPa)
205			real hybrid(size(p)) ! hybrid sigma-pressure coordinate
206
207			! Local:
208			integer l
209
210			!-------------------------------------------------------
211
212			y = p(1) / preff
213			hybrid(1) = rtsafe(funcd, x1 = 0., x2 = 1., xacc = 1e-4)
214
215			do l = 2, size(p)
216			y = p(l) / preff
217			! Assuming descending order in pressure:
218			hybrid(l) = rtsafe(funcd, x1 = 0., x2 = hybrid(l - 1), &
219			xacc = hybrid(l - 1) * 1e-4)
220			end do
221
222			end function hybrid
223
224			!**********************************************************
225
226			SUBROUTINE funcd(s, fval, fderiv)
227
228			REAL, INTENT(IN):: s
229			REAL, INTENT(OUT):: fval, fderiv
230
231			! Local:
232			real b
233
234			!------------------------------------
235
236			if (s**3 > 1. / huge(0.)) then
237			b = exp(1. - 1. / s**2)
238			fval = ya * s + (1. - ya) * b - y
239			fderiv = ya + 2 * (1. - ya) * b / s**3
240			else
241			fval = ya * s - y
242			fderiv = ya
243			end if
244
245			END SUBROUTINE funcd
246
247	guez	66	end module disvert_m