! $Id$

module o3_chem_m

  IMPLICIT none

  private o3_prod

contains

    ! This procedure evolves the ozone mass fraction through a time

    ! step taking only chemistry into account.

    ! All the 2-dimensional arrays are on the partial "physics" grid.

    ! Their shape is "(/klon, nbp_lev/)".

    ! Index "(i, :)" is for longitude "rlon(i)", latitude "rlat(i)".

    use assert_m, only: assert

    use dimphy, only: klon

    use regr_pr_comb_coefoz_m, only: c_Mob, a4_mass, a2, r_het_interm

    use mod_grid_phy_lmdz, only: nbp_lev

    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(:, :) ! (klon, nbp_lev) temperature, in K

    real, intent(in):: zmasse(:, :) ! (klon, nbp_lev)

    ! (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(in):: rlat(:), rlon(:)

    ! (longitude and latitude of each horizontal position, in degrees)

    real, intent(inout):: q(:, :) ! (klon, nbp_lev) mass fraction of ozone

    ! "q(:, k)" is at middle of layer "k".)

    ! Variables local to the procedure:

    ! (for "pi")

    integer k

    ! (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".)

    ! (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".)

    ! (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 trash1

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

    call assert(klon == (/size(q, 1), size(t_seri, 1), size(zmasse, 1), &

    call assert(nbp_lev == (/size(q, 2), size(t_seri, 2), size(zmasse, 2)/), &

    ! Compute coefficient "b":

    ! Heterogeneous chemistry is only at low temperature:

       b = r_het_interm

    elsewhere

       b = 0.

    end where

    ! Heterogeneous chemistry is only during daytime:

       where (pmu0 <= cos(87. / 180. * pi)) b(:, k) = 0.

    end forall

    ! Midpoint method:

    ! Trial step to the midpoint:

    ! "Real" step across the whole interval:

    ! Confine the mass fraction:

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

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

    ! All the 2-dimensional arrays are on the partial "physics" grid.

    ! Their shape is "(/klon, nbp_lev/)".

    ! Index "(i, :)" is for longitude "rlon(i)", latitude "rlat(i)".

    use regr_pr_comb_coefoz_m, only: a6_mass

    use assert_m, only: assert

    use dimphy, only: klon

    use mod_grid_phy_lmdz, only: nbp_lev

    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, nbp_lev)

    ! (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:

    ! (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), &

    call assert(nbp_lev == (/size(q, 2), size(zmasse, 2), size(c, 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":

    end do

  end function o3_prod

end module o3_chem_m