phylmd/Mobidic/o3_chem.f

module o3_chem_m

  ! This module is clean: no C preprocessor directive, no include line.

  IMPLICIT none

  private o3_prod

contains

  subroutine o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, q)

    ! This procedure evolves the ozone mass fraction through a time
    ! step taking only chemistry into account.

    use nrutil, only: assert
    use dimphy, only: klon
    use dimens_m, only: llm
    use regr_pr_comb_coefoz_m, only: c_Mob, a4_mass, a2, r_het_interm
    use orbite_m, only: orbite, zenang
    use nrtype, only: pi

    integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
    real, intent(in):: gmtime ! heure de la journée en fraction de jour
    real, intent(in):: t_seri(:, :) ! temperature,  in K

    real, intent(in):: zmasse(:, :)
    ! (column-density of mass of air in a cell, in kg m-2)
    ! (On the "physics" grid.
    ! "zmasse(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", for
    ! layer "k".)

    real, intent(in):: pdtphys ! time step for physics, in s

    real, intent(inout):: q(:, :) ! mass fraction of ozone
    ! (On the "physics" grid.
    ! "q(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    ! Variables local to the procedure:
    integer k

    real c(klon, llm)
    ! (constant term during a time step in the net mass production
    ! rate of ozone by chemistry, per unit mass of air, in s-1)
    ! (On the "physics" grid.
    ! "c(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    real b(klon, llm)
    ! (coefficient of "q" in the net mass production
    ! rate of ozone by chemistry, per unit mass of air, in s-1)
    ! (On the "physics" grid.
    ! "b(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    real dq_o3_chem(klon, llm)
    ! (variation of ozone mass fraction due to chemistry during a time step)
    ! (On the "physics" grid.
    ! "dq_o3_chem(i, k)" is at longitude "rlon(i)", latitude
    ! "rlat(i)", middle of layer "k".)

    real earth_long
    ! (longitude vraie de la Terre dans son orbite solaire, par
    ! rapport au point vernal (21 mars), en degrés)

    real pmu0(klon) ! mean of cosine of solar zenith angle during "pdtphys"

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

    call assert(klon == (/size(q, 1), size(t_seri, 1), size(zmasse, 1)/), &
         "o3_chem klon")
    call assert(llm == (/size(q, 2), size(t_seri, 2), size(zmasse, 2)/), &
         "o3_chem llm")

    c = c_Mob + a4_mass * t_seri

    ! Compute coefficient "b":

    ! Heterogeneous chemistry is only at low temperature:
    where (t_seri < 195.)
       b = r_het_interm
    elsewhere
       b = 0.
    end where

    ! Heterogeneous chemistry is only during daytime:
    call orbite(real(julien), earth_long)
    call zenang(earth_long, gmtime, pdtphys, pmu0)
    forall (k = 1: llm)
       where (pmu0 <= cos(87. / 180. * pi)) b(:, k) = 0.
    end forall

    b = b + a2

    ! Midpoint method:

    ! Trial step to the midpoint:
    dq_o3_chem = o3_prod(q, zmasse, c, b) * pdtphys  / 2
    ! "Real" step across the whole interval:
    dq_o3_chem = o3_prod(q + dq_o3_chem, zmasse, c, b) * pdtphys
    q = q + dq_o3_chem

    ! Confine the mass fraction:
    q = min(max(q, 0.), .01)

  end subroutine o3_chem

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

  function o3_prod(q, zmasse, c, b)

    ! This function computes the production rate of ozone by chemistry.

    use regr_pr_comb_coefoz_m, only: a6_mass
    use nrutil, only: assert
    use dimens_m, only: llm
    use dimphy, only: klon

    real, intent(in):: q(:, :) ! mass fraction of ozone
    ! (On the "physics" grid.
    ! "q(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    real, intent(in):: zmasse(:, :)
    ! (column-density of mass of air in a layer, in kg m-2)
    ! (On the "physics" grid.
    ! "zmasse(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    real, intent(in):: c(:, :)
    ! (constant term during a time step in the net mass production
    ! rate of ozone by chemistry, per unit mass of air, in s-1)
    ! (On the "physics" grid.
    ! "c(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    real, intent(in):: b(:, :)
    ! (coefficient of "q" in the net mass production
    ! rate of ozone by chemistry, per unit mass of air, in s-1)
    ! (On the "physics" grid.
    ! "b(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    real o3_prod(klon, llm)
    ! (net mass production rate of ozone by chemistry, per unit mass
    ! of air, in s-1)
    ! (On the "physics" grid.
    ! "o3_prod(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
    ! layer "k".)

    ! Variables local to the procedure:

    real sigma_mass(klon, llm)
    ! (mass column-density of ozone above point, in kg m-2)
    ! (On the "physics" grid.
    ! "sigma_mass(i, k)" is at longitude "rlon(i)", latitude
    ! "rlat(i)", middle of layer "k".)

    integer k

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

    call assert(klon == (/size(q, 1), size(zmasse, 1), size(c, 1), &
         size(b, 1)/), "o3_prod 1")
    call assert(llm == (/size(q, 2), size(zmasse, 2), size(c, 2), &
         size(b, 2)/), "o3_prod 2")

    ! Compute the column-density above the base of layer
    ! "k", and, as a first approximation, take it as column-density
    ! above the middle of layer "k":
    sigma_mass(:, llm) = zmasse(:, llm) * q(:, llm) ! top layer
    do k =  llm - 1, 1, -1
       sigma_mass(:, k) = sigma_mass(:, k+1) + zmasse(:, k) * q(:, k)
    end do

    o3_prod = c + b * q + a6_mass * sigma_mass

  end function o3_prod

end module o3_chem_m
1	guez	3	module o3_chem_m
2
3			! This module is clean: no C preprocessor directive, no include line.
4
5			IMPLICIT none
6
7			private o3_prod
8
9			contains
10
11			subroutine o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, q)
12
13			! This procedure evolves the ozone mass fraction through a time
14			! step taking only chemistry into account.
15
16			use nrutil, only: assert
17			use dimphy, only: klon
18			use dimens_m, only: llm
19	guez	9	use regr_pr_comb_coefoz_m, only: c_Mob, a4_mass, a2, r_het_interm
20	guez	3	use orbite_m, only: orbite, zenang
21			use nrtype, only: pi
22
23			integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
24			real, intent(in):: gmtime ! heure de la journée en fraction de jour
25			real, intent(in):: t_seri(:, :) ! temperature, in K
26
27			real, intent(in):: zmasse(:, :)
28			! (column-density of mass of air in a cell, in kg m-2)
29			! (On the "physics" grid.
30			! "zmasse(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", for
31			! layer "k".)
32
33			real, intent(in):: pdtphys ! time step for physics, in s
34
35			real, intent(inout):: q(:, :) ! mass fraction of ozone
36			! (On the "physics" grid.
37			! "q(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
38			! layer "k".)
39
40			! Variables local to the procedure:
41	guez	7	integer k
42	guez	3
43			real c(klon, llm)
44			! (constant term during a time step in the net mass production
45			! rate of ozone by chemistry, per unit mass of air, in s-1)
46			! (On the "physics" grid.
47			! "c(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
48			! layer "k".)
49
50			real b(klon, llm)
51			! (coefficient of "q" in the net mass production
52			! rate of ozone by chemistry, per unit mass of air, in s-1)
53			! (On the "physics" grid.
54			! "b(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
55			! layer "k".)
56
57			real dq_o3_chem(klon, llm)
58			! (variation of ozone mass fraction due to chemistry during a time step)
59			! (On the "physics" grid.
60			! "dq_o3_chem(i, k)" is at longitude "rlon(i)", latitude
61			! "rlat(i)", middle of layer "k".)
62
63			real earth_long
64			! (longitude vraie de la Terre dans son orbite solaire, par
65			! rapport au point vernal (21 mars), en degrés)
66
67			real pmu0(klon) ! mean of cosine of solar zenith angle during "pdtphys"
68
69			!-------------------------------------------------------------
70
71			call assert(klon == (/size(q, 1), size(t_seri, 1), size(zmasse, 1)/), &
72			"o3_chem klon")
73			call assert(llm == (/size(q, 2), size(t_seri, 2), size(zmasse, 2)/), &
74			"o3_chem llm")
75
76	guez	7	c = c_Mob + a4_mass * t_seri
77	guez	3
78			! Compute coefficient "b":
79
80			! Heterogeneous chemistry is only at low temperature:
81			where (t_seri < 195.)
82	guez	7	b = r_het_interm
83	guez	3	elsewhere
84			b = 0.
85			end where
86
87			! Heterogeneous chemistry is only during daytime:
88			call orbite(real(julien), earth_long)
89			call zenang(earth_long, gmtime, pdtphys, pmu0)
90			forall (k = 1: llm)
91			where (pmu0 <= cos(87. / 180. * pi)) b(:, k) = 0.
92			end forall
93
94	guez	7	b = b + a2
95	guez	3
96			! Midpoint method:
97
98			! Trial step to the midpoint:
99	guez	7	dq_o3_chem = o3_prod(q, zmasse, c, b) * pdtphys / 2
100	guez	3	! "Real" step across the whole interval:
101	guez	7	dq_o3_chem = o3_prod(q + dq_o3_chem, zmasse, c, b) * pdtphys
102	guez	3	q = q + dq_o3_chem
103
104			! Confine the mass fraction:
105			q = min(max(q, 0.), .01)
106
107			end subroutine o3_chem
108
109			!*************************************************
110
111	guez	7	function o3_prod(q, zmasse, c, b)
112	guez	3
113			! This function computes the production rate of ozone by chemistry.
114
115	guez	9	use regr_pr_comb_coefoz_m, only: a6_mass
116	guez	3	use nrutil, only: assert
117			use dimens_m, only: llm
118			use dimphy, only: klon
119
120			real, intent(in):: q(:, :) ! mass fraction of ozone
121			! (On the "physics" grid.
122			! "q(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
123			! layer "k".)
124
125			real, intent(in):: zmasse(:, :)
126			! (column-density of mass of air in a layer, in kg m-2)
127			! (On the "physics" grid.
128			! "zmasse(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
129			! layer "k".)
130
131			real, intent(in):: c(:, :)
132			! (constant term during a time step in the net mass production
133			! rate of ozone by chemistry, per unit mass of air, in s-1)
134			! (On the "physics" grid.
135			! "c(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
136			! layer "k".)
137
138			real, intent(in):: b(:, :)
139			! (coefficient of "q" in the net mass production
140			! rate of ozone by chemistry, per unit mass of air, in s-1)
141			! (On the "physics" grid.
142			! "b(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
143			! layer "k".)
144
145			real o3_prod(klon, llm)
146			! (net mass production rate of ozone by chemistry, per unit mass
147			! of air, in s-1)
148			! (On the "physics" grid.
149			! "o3_prod(i, k)" is at longitude "rlon(i)", latitude "rlat(i)", middle of
150			! layer "k".)
151
152			! Variables local to the procedure:
153
154			real sigma_mass(klon, llm)
155			! (mass column-density of ozone above point, in kg m-2)
156			! (On the "physics" grid.
157			! "sigma_mass(i, k)" is at longitude "rlon(i)", latitude
158			! "rlat(i)", middle of layer "k".)
159
160			integer k
161
162			!-------------------------------------------------------------------
163
164			call assert(klon == (/size(q, 1), size(zmasse, 1), size(c, 1), &
165			size(b, 1)/), "o3_prod 1")
166			call assert(llm == (/size(q, 2), size(zmasse, 2), size(c, 2), &
167			size(b, 2)/), "o3_prod 2")
168
169			! Compute the column-density above the base of layer
170			! "k", and, as a first approximation, take it as column-density
171			! above the middle of layer "k":
172			sigma_mass(:, llm) = zmasse(:, llm) * q(:, llm) ! top layer
173			do k = llm - 1, 1, -1
174			sigma_mass(:, k) = sigma_mass(:, k+1) + zmasse(:, k) * q(:, k)
175			end do
176
177	guez	7	o3_prod = c + b * q + a6_mass * sigma_mass
178	guez	3
179			end function o3_prod
180
181			end module o3_chem_m