source: LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/derivs_m.F90 @ 3899

module derivs_m

  implicit none

contains

  function derivs(t, y)

    use derivs_param, only: eps_cp_m, beta_m, av_thvu_constant_m, k_star_m, &
         theta_v_uns_m, theta_vs_m, rate_uns_m, av_rho_m, c_dw_m
    use nrtype, only: pi

    REAL, INTENT(IN):: t ! time, in s
    REAL, INTENT(IN):: y(:) ! (5) (/theta_v_cp, z_cp, r, av_thvu, h/)

    REAL derivs(size(y)) ! time derivative of y

    ! Local:

    real, parameter:: g = 9.81 ! acceleration of gravity, in m s-2

    real, parameter:: C_dx = 0.01
    ! drag coefficient outside cold pools, in m s-1

    real, parameter:: delta_theta = 2.
    ! jump at the top of the mixed layer, in K

    real surf_buoy_env ! surface buoyancy flux in the environment, in m s-1 K
    ! \overline{w' \theta'_{v, \mathrm{env}}}(0)

    real theta_v_cp, z_cp, r, av_thvu, h, v
    real c_star ! spreading velocity, in m s-1

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

    theta_v_cp = y(1)
    z_cp = y(2)
    r = y(3)
    av_thvu = y(4)

    v = rate_uns_m / (pi * r**2 * av_rho_m)
    c_star = k_star_m * sqrt(max(2. * g * (av_thvu - theta_v_cp) &
         / av_thvu * z_cp, 0.))

    derivs(1) = (v * (theta_v_uns_m - theta_v_cp) + c_dw_m * (theta_vs_m &
         - theta_v_cp)) / z_cp + eps_cp_m * c_star / r * (av_thvu - theta_v_cp)
    derivs(2) = v + (eps_cp_m - 2.) * z_cp / r * c_star
    derivs(3) = c_star

    if (av_thvu_constant_m) then
       derivs(4:5) = 0.
    else
       surf_buoy_env = c_dx * (theta_vs_m - av_thvu)
       h = y(5)
       derivs(4) = (1. + beta_m) / h * surf_buoy_env
       derivs(5) = beta_m / delta_theta * surf_buoy_env
    end if

  END function derivs

end module derivs_m