! $Id$
module o3_chem_m

  IMPLICIT NONE

  PRIVATE o3_prod

CONTAINS

  SUBROUTINE o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, rlat, rlon, 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 partial "physics" grid.
    ! Their shape is "(/klon, nbp_lev/)".
    ! Index "(i, :)" is for longitude "rlon(i)", latitude "rlat(i)".

    USE lmdz_assert, ONLY: assert
    USE dimphy, ONLY: klon
    USE regr_pr_comb_coefoz_m, ONLY: c_Mob, a4_mass, a2, r_het_interm
    USE lmdz_grid_phy, ONLY: nbp_lev
    USE lmdz_physical_constants, 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

    REAL c(klon, nbp_lev)
    ! (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, nbp_lev)
    ! (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, nbp_lev)
    ! (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"
    REAL trash1
    REAL trash2(klon)

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

    CALL assert(klon == (/size(q, 1), size(t_seri, 1), size(zmasse, 1), &
         size(rlat), size(rlon)/), "o3_chem klon")
    CALL assert(nbp_lev == (/size(q, 2), size(t_seri, 2), size(zmasse, 2)/), &
         "o3_chem nbp_lev")

    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, trash1)
    CALL zenang(earth_long, gmtime, 0., pdtphys, rlat, rlon, pmu0, trash2)
    forall (k = 1: nbp_lev)
       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 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 lmdz_assert, ONLY: assert
    USE dimphy, ONLY: klon
    USE lmdz_grid_phy, 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:

    REAL sigma_mass(klon, nbp_lev)
    ! (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(nbp_lev == (/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(:, nbp_lev) = zmasse(:, nbp_lev) * q(:, nbp_lev) ! top layer
    do k =  nbp_lev - 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