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	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 numer_rec, only: assert, pi
17	use dimphy, only: klon
18	use dimens_m, only: llm
19	use regr_pr_comb_coefoz_m, only: c_Mob, a4_mass, a2, r_het_interm
20	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	integer k
41
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	c = c_Mob + a4_mass * t_seri
76
77	! Compute coefficient "b":
78
79	! Heterogeneous chemistry is only at low temperature:
80	where (t_seri < 195.)
81	b = r_het_interm
82	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	b = b + a2
94
95	! Midpoint method:
96
97	! Trial step to the midpoint:
98	dq_o3_chem = o3_prod(q, zmasse, c, b) * pdtphys / 2
99	! "Real" step across the whole interval:
100	dq_o3_chem = o3_prod(q + dq_o3_chem, zmasse, c, b) * pdtphys
101	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	function o3_prod(q, zmasse, c, b)
111
112	! This function computes the production rate of ozone by chemistry.
113
114	use regr_pr_comb_coefoz_m, only: a6_mass
115	use numer_rec, only: assert
116	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	o3_prod = c + b * q + a6_mass * sigma_mass
177
178	end function o3_prod
179
180	end module o3_chem_m