Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/cold_pool_death_condition_m.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/cold_pool_death_condition_m.F90 (revision 3899) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/cold_pool_death_condition_m.F90 (revision 3899) @@ -0,0 +1,20 @@ +module cold_pool_death_condition_m + + implicit none + + real, parameter:: delta_death = - 0.1 ! in K + +contains + + logical function cold_pool_death_condition(y) + + real, intent(in):: y(:) ! (5) + ! (/theta_v_cp, z_cp, r, av_thvu, h/) + + !-------------------------------------------------- + + cold_pool_death_condition = y(1) - y(4) >= delta_death + + end function cold_pool_death_condition + +end module cold_pool_death_condition_m Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/cold_pool_death_m.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/cold_pool_death_m.F90 (revision 3899) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/cold_pool_death_m.F90 (revision 3899) @@ -0,0 +1,128 @@ +module cold_pool_death_m + + implicit none + +contains + + subroutine cold_pool_death(z_cp0, r0, h0, eps_cp, beta, av_thvu_constant, & + k_star, av_thvu, rate_uns, theta_vs, theta_v_uns, c_dw, av_rho, & + return_code, wake_lifetime, wake_final_radius, wake_final_height, t, & + yp, dydt_p, nbad) + + use cold_pool_death_condition_m, only: cold_pool_death_condition, & + delta_death + use derivs_m, only: derivs + use derivs_param, only: set_module_var + use numer_rec_95, only: rkqs, odeint_all_points + + real, intent(in):: z_cp0 ! initial height of cold pool, in m + real, intent(in):: r0 ! initial radius of cold pool, in m + real, intent(in):: h0 ! initial height of environmental boundary layer, in m + + real, intent(in):: eps_cp + ! entrainment coefficient at the border of the cold pool + + real, intent(in):: beta ! entrainment at the top of the mixed layer + logical, intent(in):: av_thvu_constant + real, intent(in):: k_star + + real, intent(in):: av_thvu ! average on cold pool height + ! of virtual potential temperature outside the cold pool, in K + + real, intent(in):: rate_uns + ! mass rate in one unsaturated downdraft, in kg s-1 + + real, intent(in):: theta_vs + ! virtual potential temperature at the surface, in K + + real, intent(in):: theta_v_uns + ! virtual potential temperature of the unsaturated downdraft, in K + + real, intent(in):: C_dw ! drag coefficient under cold pools, in m s-1 + real, intent(in):: av_rho ! average mass density in the cold pool, in kg m-3 + integer, intent(out):: return_code + REAL, intent(out):: wake_lifetime ! in s + REAL, intent(out):: wake_final_radius, wake_final_height ! in m + real, allocatable, intent(out), optional:: t(:) ! (n) time, in s + + real, allocatable, intent(out), optional:: yp(:, :) ! (5, n) + ! yp(:, i) is (/theta_v_cp, z_cp, r, av_thvu, h/) at date t(i) + + real, allocatable, intent(out), optional:: dydt_p(:, :) ! (5, n) + ! time derivative of yp, in units of yp per s + + integer, intent(out), optional:: nbad + ! number of bad (but retried and fixed) steps taken + + ! Local: + + integer n + real delta_1 ! delta theta_v at t(n - 1) + real delta_2 ! delta theta_v at t(n) + real w ! weight of interpolation between t(n - 1) and t(n) + real, parameter:: max_lifetime = 86400. * 30. ! in s + + real, allocatable:: t_local(:) ! (n) t + real, allocatable:: yp_local(:, :) ! (5, n) yp + + !---------------------------------------------------------- + + if (rate_uns < 0) then + ! No cold pool + wake_lifetime = - 4. + wake_final_radius = - 4. + wake_final_height = - 4. + return_code = 0 + else if (theta_vs < av_thvu) then + ! Cold pool temperature would not converge to environment temperature + wake_lifetime = - 3. + wake_final_radius = - 3. + wake_final_height = - 3. + return_code = 0 + else + if (theta_v_uns >= av_thvu) then + ! Bad vertical profile + wake_lifetime = - 2. + wake_final_radius = - 2. + wake_final_height = - 2. + return_code = 0 + else + call set_module_var(eps_cp, beta, av_thvu_constant, k_star, & + theta_vs, rate_uns, av_rho, theta_v_uns, c_dw) + call odeint_all_points(derivs, rkqs, t_local, yp_local, return_code, & + (/theta_v_uns, z_cp0, r0, av_thvu, h0/), nbad, dydt_p, & + x1 = 0., x2 = max_lifetime, eps = 1e-5, dmin = 0., & + stop_condition = cold_pool_death_condition) + n = size(t_local) + + if (n == 1) then + ! Dead-born + wake_lifetime = - 1. + wake_final_radius = - 1. + wake_final_height = - 1. + else if (cold_pool_death_condition(yp_local(:, n))) then + ! From the last two integration points, linearly interpolate to + ! delta theta_v = delta_death to find wake_lifetime: + delta_1 = yp_local(1, n - 1) - yp_local(4, n - 1) + delta_2 = yp_local(1, n) - yp_local(4, n) + w = (delta_death - delta_1) / (delta_2 - delta_1) + wake_lifetime = w * t_local(n) + (1. - w) * t_local(n - 1) + wake_final_radius = w * yp_local(3, n) & + + (1. - w) * yp_local(3, n - 1) + wake_final_height = w * yp_local(2, n) & + + (1. - w) * yp_local(2, n - 1) + else + ! Cold pool did not die. + wake_lifetime = - 5. + wake_final_radius = - 5. + wake_final_height = - 5. + end if + + if (present(t)) t = t_local + if (present(yp)) yp = yp_local + end if + end if + + end subroutine cold_pool_death + +end module cold_pool_death_m Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/derivs_m.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/derivs_m.F90 (revision 3899) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/derivs_m.F90 (revision 3899) @@ -0,0 +1,61 @@ +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 Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/derivs_param.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/derivs_param.F90 (revision 3899) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/Cold_pool_death/derivs_param.F90 (revision 3899) @@ -0,0 +1,72 @@ +module derivs_param + + implicit none + + real eps_cp_m ! entrainment coefficient at the border of the cold pool + !$omp THREADPRIVATE(eps_cp) + + real beta_m ! entrainment at the top of the mixed layer + !$omp THREADPRIVATE(beta) + + logical av_thvu_constant_m + !$omp THREADPRIVATE(av_thvu_constant) + + real k_star_m + !$omp THREADPRIVATE(k_star) + + real theta_vs_m ! surface virtual potential temperature, in K + !$omp THREADPRIVATE(theta_vs_m) + + real rate_uns_m ! mass rate in one unsaturated downdraft, in kg s-1 + !$omp THREADPRIVATE(rate_uns) + + real av_rho_m ! average mass density in the cold pool, in kg m-3 + !$omp THREADPRIVATE(av_rho_m) + + real theta_v_uns_m + ! virtual potential temperature of the unsaturated downdraft, in K + !$omp THREADPRIVATE(theta_v_uns_m) + + real C_dw_m ! drag coefficient under cold pools, in m s-1 + !$omp THREADPRIVATE(C_dw_m) + +contains + + subroutine set_module_var(eps_cp, beta, av_thvu_constant, k_star, theta_vs, & + rate_uns, av_rho, theta_v_uns, c_dw) + + real, intent(in):: eps_cp + ! entrainment coefficient at the border of the cold pool + + real, intent(in):: beta ! entrainment at the top of the mixed layer + logical, intent(in):: av_thvu_constant + real, intent(in):: k_star + + real, intent(in):: theta_vs + ! virtual potential temperature at the surface, in K + + real, intent(in):: rate_uns + ! mass rate in one unsaturated downdraft, in kg s-1 + + real, intent(in):: av_rho ! average mass density in the cold pool, in kg m-3 + + real, intent(in):: theta_v_uns + ! virtual potential temperature of the unsaturated downdraft, in K + + real, intent(in):: C_dw ! drag coefficient under cold pools, in m s-1 + + !---------------------------------------------- + + eps_cp_m = eps_cp + beta_m = beta + av_thvu_constant_m = av_thvu_constant + k_star_m = k_star + theta_vs_m = theta_vs + rate_uns_m = rate_uns + av_rho_m = av_rho + theta_v_uns_m = theta_v_uns + c_dw_m = c_dw + + end subroutine set_module_var + +end module derivs_param Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/calwake.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/calwake.F90 (revision 3898) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/calwake.F90 (revision 3899) @@ -14,5 +14,5 @@ wake_omg, wake_dp_deltomg, & wake_spread, wake_cstar, wake_d_deltat_gw, & - wake_ddeltat, wake_ddeltaq, wake_ds, wake_ddens) + wake_ddeltat, wake_ddeltaq, wake_ds, wake_ddens, itap) ! ************************************************************** ! * @@ -68,4 +68,5 @@ REAL, DIMENSION(klon, klev), INTENT (OUT) :: wake_ddeltat, wake_ddeltaq REAL, DIMENSION(klon), INTENT (OUT) :: wake_ds, wake_ddens + integer, intent(in):: itap ! index of call to physiq @@ -211,5 +212,5 @@ dtke, dqke, omg, dp_deltomg, spread, cstar, & d_deltat_gw, & - d_deltatw, d_deltaqw, d_sigmaw, d_wdens) ! tendencies + d_deltatw, d_deltaqw, d_sigmaw, d_wdens, itap) ! tendencies ! Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_ctrlout_mod.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 3898) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 3899) @@ -1476,4 +1476,13 @@ TYPE(ctrl_out), SAVE :: o_wake_omg = ctrl_out((/ 4, 5, 10, 10, 10, 10, 11, 11, 11, 11/), & 'wake_omg', 'wake_omg', ' ', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_wake_lifetime = ctrl_out((/ 4, 5, 10, 10, 10, 10, & + 11, 11, 11, 11/), 'wake_lifetime', 'lifetime of wake', 's', & + (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_wake_final_radius = ctrl_out((/ 4, 5, 10, 10, 10, & + 10, 11, 11, 11, 11/), 'wake_final_radius', 'final radius of wake', 'm', & + (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_wake_final_height = ctrl_out((/ 4, 5, 10, 10, 10, & + 10, 11, 11, 11, 11/), 'wake_final_height', 'final height of wake', 'm', & + (/ ('', i=1, 10) /)) TYPE(ctrl_out), SAVE :: o_wdtrainA = ctrl_out((/ 4, 5, 10, 4, 10, 10, 11, 11, 11, 11 /), & 'wdtrainA', 'precipitation from AA', '-', (/ ('', i=1, 10) /)) Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_var_mod.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_var_mod.F90 (revision 3898) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_var_mod.F90 (revision 3899) @@ -114,4 +114,10 @@ !$OMP THREADPRIVATE(sens_prec_liq_o, sens_prec_sol_o,lat_prec_liq_o,lat_prec_sol_o) + REAL, allocatable:: wake_lifetime(:) ! (klon) in s + !$OMP THREADPRIVATE(wake_lifetime) + + REAL, allocatable:: wake_final_radius(:), wake_final_height(:) ! (klon) in m + !$OMP THREADPRIVATE(wake_final_radius, wake_final_height) + CONTAINS @@ -123,4 +129,7 @@ include "clesphys.h" + + include "thermcell.h" + ! (for iflag_wake) !------------------------------------------------ @@ -163,4 +172,7 @@ IF (ok_gwd_rando) allocate(zustr_gwd_rando(klon), zvstr_gwd_rando(klon)) + if (iflag_wake>=1) allocate(wake_lifetime(klon), wake_final_radius(klon), & + wake_final_height(klon)) + END SUBROUTINE phys_output_var_init Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_write_mod.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_write_mod.F90 (revision 3898) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/phys_output_write_mod.F90 (revision 3899) @@ -182,5 +182,5 @@ o_map_emis_Anv, o_map_pcld_Anv, o_map_tcld_Anv, & o_map_ntot, o_map_hc,o_map_hist,o_map_Cb,o_map_ThCi,o_map_Anv, & - o_alt_tropo + o_alt_tropo, o_wake_lifetime, o_wake_final_radius, o_wake_final_height #ifdef CPP_StratAer @@ -317,5 +317,5 @@ map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,& map_emis_Anv,map_pcld_Anv,map_tcld_Anv, & - alt_tropo + alt_tropo, wake_lifetime, wake_final_radius, wake_final_height USE ocean_slab_mod, ONLY: nslay, tslab, slab_bils, slab_bilg, tice, & @@ -1057,4 +1057,7 @@ CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d) CALL histwrite_phy(o_dqwak2d, zx_tmp_fi2d) + CALL histwrite_phy(o_wake_lifetime, wake_lifetime) + CALL histwrite_phy(o_wake_final_radius, wake_final_radius) + CALL histwrite_phy(o_wake_final_height, wake_final_height) ENDIF ! iflag_wake>=1 CALL histwrite_phy(o_ftd, ftd) Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/physiq_mod.F90 (revision 3898) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/physiq_mod.F90 (revision 3899) @@ -2772,5 +2772,5 @@ wake_omg, wake_dp_deltomg, & wake_spread, wake_Cstar, d_deltat_wk_gw, & - d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk) + d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, itap) ! !jyg Reinitialize itapwk when wakes have been called Index: LMDZ5/branches/Cold_pool_death/libf/phylmd/wake.F90 =================================================================== --- LMDZ5/branches/Cold_pool_death/libf/phylmd/wake.F90 (revision 3898) +++ LMDZ5/branches/Cold_pool_death/libf/phylmd/wake.F90 (revision 3899) @@ -11,5 +11,5 @@ dtke, dqke, omg, dp_deltomg, spread, cstar, & d_deltat_gw, & - d_deltatw2, d_deltaqw2, d_sigmaw2, d_wdens2) ! tendencies + d_deltatw2, d_deltaqw2, d_sigmaw2, d_wdens2, itap) ! tendencies @@ -26,5 +26,11 @@ USE dimphy use mod_phys_lmdz_para + use cold_pool_death_m, only: cold_pool_death USE print_control_mod, ONLY: prt_level + use write_field_phy, only: WriteField_phy + use phys_local_var_mod, only: zxtsol, q_seri, sens_w + use phys_output_var_mod, only: wake_lifetime, wake_final_radius, & + wake_final_height + IMPLICIT NONE ! ============================================================================ @@ -158,4 +164,5 @@ REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2 REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_wdens2 + integer, intent(in):: itap ! index of call to physiq ! Variables internes @@ -272,4 +279,33 @@ REAL, DIMENSION(klon) :: sigmaw_in ! pour les prints + + real:: z_cp0 = 1000. ! initial height of cold pool, in m + !$omp THREADPRIVATE(z_cp0) + + real:: r0 = 2000. ! initial radius of cold pool, in m + !$omp THREADPRIVATE(r0) + + real:: h0 = 1000. ! initial height of environmental boundary layer, in m + !$omp THREADPRIVATE(h0) + + real:: eps_cp = 0.06 ! entrainment coefficient at the border of the cold pool + !$omp THREADPRIVATE(eps_cp) + + real:: beta = 0.2 ! entrainment at the top of the mixed layer + !$omp THREADPRIVATE(beta) + + logical:: av_thvu_constant = .true. + !$omp THREADPRIVATE(av_thvu_constant) + + real:: k_star = 0.68 + !$omp THREADPRIVATE(k_star) + + real rate_uns(klon) ! mass rate in one unsaturated downdraft, in kg s-1 + real theta_vs(klon) ! virtual potential temperature at the surface, in K + + real theta_v_uns(klon) + ! virtual potential temperature of the unsaturated downdraft, in K + + integer return_code(klon) ! ------------------------------------------------------------------------- @@ -348,4 +384,11 @@ WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl + CALL getin_p('z_cp0', z_cp0) + CALL getin_p('r0', r0) + CALL getin_p('h0', h0) + CALL getin_p('eps_cp', eps_cp) + CALL getin_p('beta', beta) + CALL getin_p('av_thvu_constant', av_thvu_constant) + CALL getin_p('k_star', k_star) first=.false. endif @@ -1970,5 +2013,37 @@ !>jyg - + do i = 1, klon + if (ktopw(i) /= 0) then + rate_uns(i) = - amdwn(i, ktopw(i)) / wdens(i) + else + rate_uns(i) = -1. + end if + end do + + theta_vs = zxtsol * (1e5 / ph(:, 1))**rkappa * (1 + epsim1 * q_seri(:, 1)) + + forall (i = 1:klon, ktopw(i) /= 0 .and. theta_vs(i) >= av_thvu(i)) & + theta_v_uns(i) = the(i, ktopw(i)) * (1. + epsim1 * qe(i, ktopw(i))) + + call WriteField_phy("thu", thu, klev) + call WriteField_phy("qu", qu, klev) + call WriteField_phy("av_thvu", av_thvu, 1) + call WriteField_phy("ktopw", real(ktopw), 1) + call WriteField_phy("wdens", wdens, 1) + call WriteField_phy("av_rho", av_rho, 1) + call WriteField_phy("rate_uns", rate_uns, 1) + call WriteField_phy("theta_vs", theta_vs, 1) + call WriteField_phy("the", the, klev) + call WriteField_phy("qe", qe, klev) + + do i = 1, klon + call cold_pool_death(z_cp0, r0, h0, eps_cp, beta, av_thvu_constant, & + k_star, av_thvu(i), rate_uns(i), theta_vs(i), theta_v_uns(i), & + sens_w(i) / (av_rho(i) * rcpd * (theta_vs(i) - theta_v_uns(i))), & + av_rho(i), return_code(i), wake_lifetime(i), wake_final_radius(i), & + wake_final_height(i)) + end do + if (any(return_code /= 0)) print *, & + "cold_pool_death: return_code /= 0 somewhere, itap = ", itap RETURN