--- trunk/Sources/phylmd/CV30_routines/cv30_yield.f 2016/03/29 15:20:23 189 +++ trunk/Sources/phylmd/CV30_routines/cv30_yield.f 2016/06/06 17:42:15 201 @@ -5,29 +5,46 @@ contains SUBROUTINE cv30_yield(icb, inb, delt, t, rr, u, v, gz, p, ph, h, hp, lv, & - cpn, th, ep, clw, m, tp, mp, rp, up, vp, wt, water, evap, b, ment, & + cpn, th, ep, clw, m, tp, mp, qp, up, vp, wt, water, evap, b, ment, & qent, uent, vent, nent, elij, sig, tv, tvp, iflag, precip, VPrecip, & ft, fr, fu, fv, upwd, dnwd, dnwd0, ma, mike, tls, tps, qcondc) + ! Tendencies, precipitation, variables of interface with other + ! processes, etc. + use conema3_m, only: iflag_clw - use cv30_param_m, only: delta, minorig, nl, sigd - use cvthermo, only: cl, cpd, cpv, grav, rowl, rrd, rrv + use cv30_param_m, only: minorig, nl, sigd + use cv_thermo_m, only: rowl USE dimphy, ONLY: klev, klon + use SUPHEC_M, only: rg, rcpd, rcw, rcpv, rd, rv ! inputs: - integer, intent(in):: icb(:), inb(:) ! (ncum) + + integer, intent(in):: icb(:) + + integer, intent(in):: inb(:) ! (ncum) + ! first model level above the level of neutral buoyancy of the + ! parcel (1 <= inb <= nl - 1) + real, intent(in):: delt - real t(klon, klev), rr(klon, klev), u(klon, klev), v(klon, klev) + real, intent(in):: t(klon, klev), rr(klon, klev) + real, intent(in):: u(klon, klev), v(klon, klev) real gz(klon, klev) real p(klon, klev) real ph(klon, klev + 1), h(klon, klev), hp(klon, klev) - real lv(klon, klev), cpn(klon, klev) - real th(klon, klev) + real, intent(in):: lv(:, :) ! (klon, klev) + + real, intent(in):: cpn(:, :) ! (ncum, nl) + ! specific heat capacity at constant pressure of humid air, in J K-1 kg-1 + + real, intent(in):: th(:, :) ! (ncum, nl) real ep(klon, klev), clw(klon, klev) real m(klon, klev) real tp(klon, klev) - real mp(klon, klev), rp(klon, klev), up(klon, klev) - real vp(klon, klev), wt(klon, klev) + real, intent(in):: mp(:, :) ! (ncum, nl) + real, intent(in):: qp(:, :), up(:, :) ! (klon, klev) + real, intent(in):: vp(:, 2:) ! (ncum, 2:nl) + real, intent(in):: wt(:, :) ! (ncum, nl - 1) real, intent(in):: water(:, :), evap(:, :) ! (ncum, nl) real, intent(in):: b(:, :) ! (ncum, nl - 1) real ment(klon, klev, klev), qent(klon, klev, klev), uent(klon, klev, klev) @@ -37,10 +54,8 @@ real sig(klon, klev) real tv(klon, klev), tvp(klon, klev) - ! input / output: - integer iflag(klon) - ! outputs: + integer, intent(out):: iflag(:) ! (ncum) real precip(klon) real VPrecip(klon, klev + 1) real ft(klon, klev), fr(klon, klev), fu(klon, klev), fv(klon, klev) @@ -52,9 +67,10 @@ real qcondc(klon, klev) ! Local: + real, parameter:: delta = 0.01 ! interface cloud parameterization integer ncum - integer i, k, il, n, j, num1 - real rat, awat, delti + integer i, k, il, n, j + real awat, delti real ax, bx, cx, dx real cpinv, rdcp, dpinv real lvcp(klon, klev) @@ -67,6 +83,7 @@ !------------------------------------------------------------- ncum = size(icb) + iflag = 0 ! initialization: @@ -96,11 +113,11 @@ enddo enddo - ! calculate surface precipitation in mm / day + ! calculate surface precipitation in mm / day do il = 1, ncum if (ep(il, inb(il)) >= 1e-4) precip(il) = wt(il, 1) * sigd & - * water(il, 1) * 86400. * 1000. / (rowl * grav) + * water(il, 1) * 86400. * 1000. / (rowl * rg) enddo ! CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg / m2 / s === @@ -109,12 +126,12 @@ do k = 1, nl - 1 do il = 1, ncum if (k <= inb(il)) VPrecip(il, k) = wt(il, k) * sigd * water(il, k) & - / grav + / rg end do end do - ! calculate tendencies of lowest level potential temperature - ! and mixing ratio + ! calculate tendencies of lowest level potential temperature + ! and mixing ratio do il = 1, ncum work(il) = 1.0 / (ph(il, 1) - ph(il, 2)) @@ -128,63 +145,51 @@ enddo do il = 1, ncum - ! Consist vect: - if (0.01 * grav * work(il) * am(il) >= delti) iflag(il) = 1 - - ft(il, 1) = 0.01 * grav * work(il) * am(il) * (t(il, 2) - t(il, 1) & - + (gz(il, 2) - gz(il, 1)) / cpn(il, 1)) + if (0.01 * rg * work(il) * am(il) >= delti) iflag(il) = 1 - ft(il, 1) = ft(il, 1) - 0.5 * lvcp(il, 1) * sigd * (evap(il, 1) & - + evap(il, 2)) - - ft(il, 1) = ft(il, 1) - 0.009 * grav * sigd * mp(il, 2) & - * t(il, 1) * b(il, 1) * work(il) - - ft(il, 1) = ft(il, 1) + 0.01 * sigd * wt(il, 1) * (cl - cpd) & - * water(il, 2) * (t(il, 2) - t(il, 1)) * work(il) / cpn(il, 1) + ft(il, 1) = 0.01 * rg * work(il) * am(il) * (t(il, 2) - t(il, 1) & + + (gz(il, 2) - gz(il, 1)) / cpn(il, 1)) - 0.5 * lvcp(il, 1) & + * sigd * (evap(il, 1) + evap(il, 2)) - 0.009 * rg * sigd & + * mp(il, 2) * t(il, 1) * b(il, 1) * work(il) + 0.01 * sigd & + * wt(il, 1) * (rcw - rcpd) * water(il, 2) * (t(il, 2) - t(il, 1)) & + * work(il) / cpn(il, 1) !jyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3) - ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap) - fr(il, 1) = 0.01 * grav * mp(il, 2) * (rp(il, 2) - rr(il, 1)) & + ! (sb: pour l'instant, on ne fait que le chgt concernant rg, pas evap) + fr(il, 1) = 0.01 * rg * mp(il, 2) * (qp(il, 2) - rr(il, 1)) & * work(il) + sigd * 0.5 * (evap(il, 1) + evap(il, 2)) ! + tard : + sigd * evap(il, 1) - fr(il, 1) = fr(il, 1) + 0.01 * grav * am(il) * (rr(il, 2) - rr(il, 1)) & + fr(il, 1) = fr(il, 1) + 0.01 * rg * am(il) * (rr(il, 2) - rr(il, 1)) & * work(il) - fu(il, 1) = fu(il, 1) + 0.01 * grav * work(il) * (mp(il, 2) & + fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) & * (up(il, 2) - u(il, 1)) + am(il) * (u(il, 2) - u(il, 1))) - fv(il, 1) = fv(il, 1) + 0.01 * grav * work(il) * (mp(il, 2) & + fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * (mp(il, 2) & * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1))) - enddo ! il + enddo do j = 2, nl do il = 1, ncum if (j <= inb(il)) then - fr(il, 1) = fr(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) & + fr(il, 1) = fr(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) & * (qent(il, j, 1) - rr(il, 1)) - fu(il, 1) = fu(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) & + fu(il, 1) = fu(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) & * (uent(il, j, 1) - u(il, 1)) - fv(il, 1) = fv(il, 1) + 0.01 * grav * work(il) * ment(il, j, 1) & + fv(il, 1) = fv(il, 1) + 0.01 * rg * work(il) * ment(il, j, 1) & * (vent(il, j, 1) - v(il, 1)) endif enddo enddo - ! calculate tendencies of potential temperature and mixing ratio - ! at levels above the lowest level + ! calculate tendencies of potential temperature and mixing ratio + ! at levels above the lowest level - ! first find the net saturated updraft and downdraft mass fluxes - ! through each level + ! first find the net saturated updraft and downdraft mass fluxes + ! through each level loop_i: do i = 2, nl - 1 - num1 = 0 - - do il = 1, ncum - if (i <= inb(il)) num1 = num1 + 1 - enddo - - if (num1 > 0) then + if (any(inb >= i)) then amp1(:ncum) = 0. ad(:ncum) = 0. @@ -221,32 +226,26 @@ dpinv = 1.0 / (ph(il, i) - ph(il, i + 1)) cpinv = 1.0 / cpn(il, i) - ! Vecto: - if (0.01 * grav * dpinv * amp1(il) >= delti) iflag(il) = 1 + if (0.01 * rg * dpinv * amp1(il) >= delti) iflag(il) = 1 - ft(il, i) = 0.01 * grav * dpinv * (amp1(il) * (t(il, i + 1) & + ft(il, i) = 0.01 * rg * dpinv * (amp1(il) * (t(il, i + 1) & - t(il, i) + (gz(il, i + 1) - gz(il, i)) * cpinv) & - ad(il) * (t(il, i) - t(il, i - 1) + (gz(il, i) & - gz(il, i - 1)) * cpinv)) - 0.5 * sigd * lvcp(il, i) & - * (evap(il, i) + evap(il, i + 1)) - rat = cpn(il, i - 1) * cpinv - ft(il, i) = ft(il, i) - 0.009 * grav * sigd * (mp(il, i + 1) & - * t(il, i) * b(il, i) - mp(il, i) * t(il, i - 1) * rat & - * b(il, i - 1)) * dpinv - ft(il, i) = ft(il, i) + 0.01 * grav * dpinv * ment(il, i, i) & - * (hp(il, i) - h(il, i) + t(il, i) * (cpv - cpd) & - * (rr(il, i) - qent(il, i, i))) * cpinv - - ft(il, i) = ft(il, i) + 0.01 * sigd * wt(il, i) * (cl - cpd) & - * water(il, i + 1) * (t(il, i + 1) - t(il, i)) * dpinv & - * cpinv - - fr(il, i) = 0.01 * grav * dpinv * (amp1(il) * (rr(il, i + 1) & + * (evap(il, i) + evap(il, i + 1)) - 0.009 * rg * sigd & + * (mp(il, i + 1) * t(il, i) * b(il, i) - mp(il, i) & + * t(il, i - 1) * cpn(il, i - 1) * cpinv * b(il, i - 1)) & + * dpinv + 0.01 * rg * dpinv * ment(il, i, i) & + * (hp(il, i) - h(il, i) + t(il, i) * (rcpv - rcpd) & + * (rr(il, i) - qent(il, i, i))) * cpinv + 0.01 * sigd & + * wt(il, i) * (rcw - rcpd) * water(il, i + 1) & + * (t(il, i + 1) - t(il, i)) * dpinv * cpinv + fr(il, i) = 0.01 * rg * dpinv * (amp1(il) * (rr(il, i + 1) & - rr(il, i)) - ad(il) * (rr(il, i) - rr(il, i - 1))) - fu(il, i) = fu(il, i) + 0.01 * grav * dpinv * (amp1(il) & + fu(il, i) = fu(il, i) + 0.01 * rg * dpinv * (amp1(il) & * (u(il, i + 1) - u(il, i)) - ad(il) * (u(il, i) & - u(il, i - 1))) - fv(il, i) = fv(il, i) + 0.01 * grav * dpinv * (amp1(il) & + fv(il, i) = fv(il, i) + 0.01 * rg * dpinv * (amp1(il) & * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) & - v(il, i - 1))) endif @@ -261,11 +260,11 @@ awat = elij(il, k, i) - (1. - ep(il, i)) * clw(il, i) awat = amax1(awat, 0.0) - fr(il, i) = fr(il, i) + 0.01 * grav * dpinv & + fr(il, i) = fr(il, i) + 0.01 * rg * dpinv & * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i)) - fu(il, i) = fu(il, i) + 0.01 * grav * dpinv & + fu(il, i) = fu(il, i) + 0.01 * rg * dpinv & * ment(il, k, i) * (uent(il, k, i) - u(il, i)) - fv(il, i) = fv(il, i) + 0.01 * grav * dpinv & + fv(il, i) = fv(il, i) + 0.01 * rg * dpinv & * ment(il, k, i) * (vent(il, k, i) - v(il, i)) ! (saturated updrafts resulting from mixing) @@ -281,11 +280,11 @@ dpinv = 1.0 / (ph(il, i) - ph(il, i + 1)) cpinv = 1.0 / cpn(il, i) - fr(il, i) = fr(il, i) + 0.01 * grav * dpinv & + fr(il, i) = fr(il, i) + 0.01 * rg * dpinv & * ment(il, k, i) * (qent(il, k, i) - rr(il, i)) - fu(il, i) = fu(il, i) + 0.01 * grav * dpinv & + fu(il, i) = fu(il, i) + 0.01 * rg * dpinv & * ment(il, k, i) * (uent(il, k, i) - u(il, i)) - fv(il, i) = fv(il, i) + 0.01 * grav * dpinv & + fv(il, i) = fv(il, i) + 0.01 * rg * dpinv & * ment(il, k, i) * (vent(il, k, i) - v(il, i)) endif end do @@ -299,14 +298,14 @@ ! sb: on ne fait pas encore la correction permettant de mieux ! conserver l'eau: fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(il, i) & - + evap(il, i + 1)) + 0.01 * grav * (mp(il, i + 1) & - * (rp(il, i + 1) - rr(il, i)) - mp(il, i) * (rp(il, i) & + + evap(il, i + 1)) + 0.01 * rg * (mp(il, i + 1) & + * (qp(il, i + 1) - rr(il, i)) - mp(il, i) * (qp(il, i) & - rr(il, i - 1))) * dpinv - fu(il, i) = fu(il, i) + 0.01 * grav * (mp(il, i + 1) & + fu(il, i) = fu(il, i) + 0.01 * rg * (mp(il, i + 1) & * (up(il, i + 1) - u(il, i)) - mp(il, i) * (up(il, i) & - u(il, i - 1))) * dpinv - fv(il, i) = fv(il, i) + 0.01 * grav * (mp(il, i + 1) & + fv(il, i) = fv(il, i) + 0.01 * rg * (mp(il, i + 1) & * (vp(il, i + 1) - v(il, i)) - mp(il, i) * (vp(il, i) & - v(il, i - 1))) * dpinv endif @@ -340,13 +339,13 @@ end if end do loop_i - ! move the detrainment at level inb down to level inb - 1 - ! in such a way as to preserve the vertically - ! integrated enthalpy and water tendencies + ! move the detrainment at level inb down to level inb - 1 + ! in such a way as to preserve the vertically + ! integrated enthalpy and water tendencies do il = 1, ncum ax = 0.1 * ment(il, inb(il), inb(il)) * (hp(il, inb(il)) & - - h(il, inb(il)) + t(il, inb(il)) * (cpv - cpd) & + - h(il, inb(il)) + t(il, inb(il)) * (rcpv - rcpd) & * (rr(il, inb(il)) - qent(il, inb(il), inb(il)))) & / (cpn(il, inb(il)) * (ph(il, inb(il)) - ph(il, inb(il) + 1))) ft(il, inb(il)) = ft(il, inb(il)) - ax @@ -377,7 +376,7 @@ end do - ! homoginize tendencies below cloud base + ! homoginize tendencies below cloud base do il = 1, ncum asum(il) = 0.0 @@ -390,9 +389,9 @@ do il = 1, ncum if (i <= (icb(il) - 1)) then asum(il) = asum(il) + ft(il, i) * (ph(il, i) - ph(il, i + 1)) - bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (cl - cpd) & + bsum(il) = bsum(il) + fr(il, i) * (lv(il, i) + (rcw - rcpd) & * (t(il, i) - t(il, 1))) * (ph(il, i) - ph(il, i + 1)) - csum(il) = csum(il) + (lv(il, i) + (cl - cpd) * (t(il, i) & + csum(il) = csum(il) + (lv(il, i) + (rcw - rcpd) * (t(il, i) & - t(il, 1))) * (ph(il, i) - ph(il, i + 1)) dsum(il) = dsum(il) + t(il, i) * (ph(il, i) - ph(il, i + 1)) & / th(il, i) @@ -409,7 +408,7 @@ enddo enddo - ! reset counter and return + ! reset counter and return do il = 1, ncum sig(il, klev) = 2.0 @@ -475,10 +474,8 @@ enddo enddo - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! determination de la variation de flux ascendant entre - ! deux niveau non dilue mike - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! D\'etermination de la variation de flux ascendant entre + ! deux niveaux non dilu\'es mike do i = 1, nl do il = 1, ncum @@ -520,23 +517,19 @@ enddo enddo - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! icb represente de niveau ou se trouve la - ! base du nuage, et inb le top du nuage - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + ! icb repr\'esente le niveau o\`u se trouve la base du nuage, et + ! inb le sommet du nuage do i = 1, klev DO il = 1, ncum - rdcp = (rrd * (1. - rr(il, i)) - rr(il, i) * rrv) & - / (cpd * (1. - rr(il, i)) + rr(il, i) * cpv) + rdcp = (rd * (1. - rr(il, i)) - rr(il, i) * rv) & + / (rcpd * (1. - rr(il, i)) + rr(il, i) * rcpv) tls(il, i) = t(il, i) * (1000.0 / p(il, i))**rdcp tps(il, i) = tp(il, i) end DO enddo - ! diagnose the in-cloud mixing ratio - ! of condensed water - ! + ! Diagnose the in-cloud mixing ratio of condensed water do i = 1, klev do il = 1, ncum @@ -562,7 +555,7 @@ do il = 1, ncum if (i >= icb(il) .and. i <= (inb(il) - 1) & .and. j >= icb(il)) then - sax(il, i) = sax(il, i) + rrd * (tvp(il, j) - tv(il, j)) & + sax(il, i) = sax(il, i) + rd * (tvp(il, j) - tv(il, j)) & * (ph(il, j) - ph(il, j + 1)) / p(il, j) endif enddo @@ -580,7 +573,7 @@ do i = 1, nl do il = 1, ncum - if (wa(il, i) > 0.0) siga(il, i) = mac(il, i) / wa(il, i) * rrd & + if (wa(il, i) > 0.0) siga(il, i) = mac(il, i) / wa(il, i) * rd & * tvp(il, i) / p(il, i) / 100. / delta siga(il, i) = min(siga(il, i), 1.0)