phylmd/Mobidic/o3_chem.f

module o3_chem_m

  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.

    ! All the 2-dimensional arrays are on the "physics" grid.
    ! Their shape is "(/klon, llm/)".
    ! Index "(i, :)" is for longitude "rlon(i)", latitude "rlat(i)".

    use nr_util, only: assert, pi
    use dimphy, only: klon
    use dimensions, only: llm
    use regr_pr_comb_coefoz_m, only: c_Mob, a4_mass, a2, r_het_interm
    use orbite_m, only: orbite
    use zenang_m, only: 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(:, :) ! (klon, llm) temperature, in K

    real, intent(in):: zmasse(:, :) ! (klon, llm)
    ! (column-density of mass of air in a cell, in kg m-2)
    ! "zmasse(:, k)" is for layer "k".)

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

    real, intent(inout):: q(:, :) ! (klon, llm) mass fraction of ozone
    ! "q(:, k)" is at 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)
    ! "c(:, k)" is at 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)
    ! "b(:, k)" is at middle of layer "k".)

    real dq_o3_chem(klon, llm)
    ! (variation of ozone mass fraction due to chemistry during a time step)
    ! "dq_o3_chem(:, k)" is at 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.

    ! All the 2-dimensional arrays are on the "physics" grid.
    ! Their shape is "(/klon, llm/)".
    ! Index "(i, :)" is for longitude "rlon(i)", latitude "rlat(i)".

    use regr_pr_comb_coefoz_m, only: a6_mass
    use nr_util, only: assert
    use dimensions, only: llm
    use dimphy, only: klon

    real, intent(in):: q(:, :) ! mass fraction of ozone
    ! "q(:, k)" is at middle of layer "k".)

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