phylmd/Mobidic/o3_chem.f90

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 numer_rec, only: assert, pi
    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

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