| 1 | module cold_pool_death_m |
|---|
| 2 | |
|---|
| 3 | implicit none |
|---|
| 4 | |
|---|
| 5 | contains |
|---|
| 6 | |
|---|
| 7 | subroutine cold_pool_death(z_cp0, r0, h0, eps_cp, beta, av_thvu_constant, & |
|---|
| 8 | k_star, av_thvu, rate_uns, theta_vs, theta_v_uns, c_dw, av_rho, & |
|---|
| 9 | return_code, wake_lifetime, wake_final_radius, wake_final_height, t, & |
|---|
| 10 | yp, dydt_p, nbad) |
|---|
| 11 | |
|---|
| 12 | use cold_pool_death_condition_m, only: cold_pool_death_condition, & |
|---|
| 13 | delta_death |
|---|
| 14 | use derivs_m, only: derivs |
|---|
| 15 | use derivs_param, only: set_module_var |
|---|
| 16 | use numer_rec_95, only: rkqs, odeint_all_points |
|---|
| 17 | |
|---|
| 18 | real, intent(in):: z_cp0 ! initial height of cold pool, in m |
|---|
| 19 | real, intent(in):: r0 ! initial radius of cold pool, in m |
|---|
| 20 | real, intent(in):: h0 ! initial height of environmental boundary layer, in m |
|---|
| 21 | |
|---|
| 22 | real, intent(in):: eps_cp |
|---|
| 23 | ! entrainment coefficient at the border of the cold pool |
|---|
| 24 | |
|---|
| 25 | real, intent(in):: beta ! entrainment at the top of the mixed layer |
|---|
| 26 | logical, intent(in):: av_thvu_constant |
|---|
| 27 | real, intent(in):: k_star |
|---|
| 28 | |
|---|
| 29 | real, intent(in):: av_thvu ! average on cold pool height |
|---|
| 30 | ! of virtual potential temperature outside the cold pool, in K |
|---|
| 31 | |
|---|
| 32 | real, intent(in):: rate_uns |
|---|
| 33 | ! mass rate in one unsaturated downdraft, in kg s-1 |
|---|
| 34 | |
|---|
| 35 | real, intent(in):: theta_vs |
|---|
| 36 | ! virtual potential temperature at the surface, in K |
|---|
| 37 | |
|---|
| 38 | real, intent(in):: theta_v_uns |
|---|
| 39 | ! virtual potential temperature of the unsaturated downdraft, in K |
|---|
| 40 | |
|---|
| 41 | real, intent(in):: C_dw ! drag coefficient under cold pools, in m s-1 |
|---|
| 42 | real, intent(in):: av_rho ! average mass density in the cold pool, in kg m-3 |
|---|
| 43 | integer, intent(out):: return_code |
|---|
| 44 | REAL, intent(out):: wake_lifetime ! in s |
|---|
| 45 | REAL, intent(out):: wake_final_radius, wake_final_height ! in m |
|---|
| 46 | real, allocatable, intent(out), optional:: t(:) ! (n) time, in s |
|---|
| 47 | |
|---|
| 48 | real, allocatable, intent(out), optional:: yp(:, :) ! (5, n) |
|---|
| 49 | ! yp(:, i) is (/theta_v_cp, z_cp, r, av_thvu, h/) at date t(i) |
|---|
| 50 | |
|---|
| 51 | real, allocatable, intent(out), optional:: dydt_p(:, :) ! (5, n) |
|---|
| 52 | ! time derivative of yp, in units of yp per s |
|---|
| 53 | |
|---|
| 54 | integer, intent(out), optional:: nbad |
|---|
| 55 | ! number of bad (but retried and fixed) steps taken |
|---|
| 56 | |
|---|
| 57 | ! Local: |
|---|
| 58 | |
|---|
| 59 | integer n |
|---|
| 60 | real delta_1 ! delta theta_v at t(n - 1) |
|---|
| 61 | real delta_2 ! delta theta_v at t(n) |
|---|
| 62 | real w ! weight of interpolation between t(n - 1) and t(n) |
|---|
| 63 | real, parameter:: max_lifetime = 86400. * 30. ! in s |
|---|
| 64 | |
|---|
| 65 | real, allocatable:: t_local(:) ! (n) t |
|---|
| 66 | real, allocatable:: yp_local(:, :) ! (5, n) yp |
|---|
| 67 | |
|---|
| 68 | !---------------------------------------------------------- |
|---|
| 69 | |
|---|
| 70 | if (rate_uns < 0) then |
|---|
| 71 | ! No cold pool |
|---|
| 72 | wake_lifetime = - 4. |
|---|
| 73 | wake_final_radius = - 4. |
|---|
| 74 | wake_final_height = - 4. |
|---|
| 75 | return_code = 0 |
|---|
| 76 | else if (theta_vs < av_thvu) then |
|---|
| 77 | ! Cold pool temperature would not converge to environment temperature |
|---|
| 78 | wake_lifetime = - 3. |
|---|
| 79 | wake_final_radius = - 3. |
|---|
| 80 | wake_final_height = - 3. |
|---|
| 81 | return_code = 0 |
|---|
| 82 | else |
|---|
| 83 | if (theta_v_uns >= av_thvu) then |
|---|
| 84 | ! Bad vertical profile |
|---|
| 85 | wake_lifetime = - 2. |
|---|
| 86 | wake_final_radius = - 2. |
|---|
| 87 | wake_final_height = - 2. |
|---|
| 88 | return_code = 0 |
|---|
| 89 | else |
|---|
| 90 | call set_module_var(eps_cp, beta, av_thvu_constant, k_star, & |
|---|
| 91 | theta_vs, rate_uns, av_rho, theta_v_uns, c_dw) |
|---|
| 92 | call odeint_all_points(derivs, rkqs, t_local, yp_local, return_code, & |
|---|
| 93 | (/theta_v_uns, z_cp0, r0, av_thvu, h0/), nbad, dydt_p, & |
|---|
| 94 | x1 = 0., x2 = max_lifetime, eps = 1e-5, dmin = 0., & |
|---|
| 95 | stop_condition = cold_pool_death_condition) |
|---|
| 96 | n = size(t_local) |
|---|
| 97 | |
|---|
| 98 | if (n == 1) then |
|---|
| 99 | ! Dead-born |
|---|
| 100 | wake_lifetime = - 1. |
|---|
| 101 | wake_final_radius = - 1. |
|---|
| 102 | wake_final_height = - 1. |
|---|
| 103 | else if (cold_pool_death_condition(yp_local(:, n))) then |
|---|
| 104 | ! From the last two integration points, linearly interpolate to |
|---|
| 105 | ! delta theta_v = delta_death to find wake_lifetime: |
|---|
| 106 | delta_1 = yp_local(1, n - 1) - yp_local(4, n - 1) |
|---|
| 107 | delta_2 = yp_local(1, n) - yp_local(4, n) |
|---|
| 108 | w = (delta_death - delta_1) / (delta_2 - delta_1) |
|---|
| 109 | wake_lifetime = w * t_local(n) + (1. - w) * t_local(n - 1) |
|---|
| 110 | wake_final_radius = w * yp_local(3, n) & |
|---|
| 111 | + (1. - w) * yp_local(3, n - 1) |
|---|
| 112 | wake_final_height = w * yp_local(2, n) & |
|---|
| 113 | + (1. - w) * yp_local(2, n - 1) |
|---|
| 114 | else |
|---|
| 115 | ! Cold pool did not die. |
|---|
| 116 | wake_lifetime = - 5. |
|---|
| 117 | wake_final_radius = - 5. |
|---|
| 118 | wake_final_height = - 5. |
|---|
| 119 | end if |
|---|
| 120 | |
|---|
| 121 | if (present(t)) t = t_local |
|---|
| 122 | if (present(yp)) yp = yp_local |
|---|
| 123 | end if |
|---|
| 124 | end if |
|---|
| 125 | |
|---|
| 126 | end subroutine cold_pool_death |
|---|
| 127 | |
|---|
| 128 | end module cold_pool_death_m |
|---|
