--- trunk/Sources/phylmd/CV30_routines/cv30_undilute2.f 2016/05/12 13:00:07 192 +++ trunk/phylmd/CV30_routines/cv30_undilute2.f 2018/02/05 10:39:38 254 @@ -4,9 +4,8 @@ contains - SUBROUTINE cv30_undilute2(ncum, icb, icbs, nk, tnk, qnk, gznk, t, qs, gz, & - p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, ep, sigp, & - buoy) + SUBROUTINE cv30_undilute2(icb, icbs, tnk, qnk, gznk, t, qs, gz, p, h, tv, & + lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, ep, buoy) ! Undilute (adiabatic) updraft, second part. Purpose: find the ! rest of the lifted parcel temperatures; compute the @@ -15,22 +14,21 @@ ! Vertical profile of buoyancy computed here (use of buoybase). - use conema3_m, only: epmax - use cv30_param_m, only: dtovsh, minorig, nl, pbcrit, ptcrit, spfac - use cv_thermo_m, only: cl, clmcpv, cpd, cpv, eps, lv0, rrv + use conf_phys_m, only: epmax + use cv30_param_m, only: minorig, nl + use cv_thermo_m, only: clmcpv, eps USE dimphy, ONLY: klon, klev + use SUPHEC_M, only: rcw, rlvtt, rcpd, rcpv, rv - integer, intent(in):: ncum - - integer, intent(in):: icb(klon), icbs(klon) + integer, intent(in):: icb(:), icbs(:) ! (ncum) ! icbs is the first level above LCL (may differ from icb) - integer, intent(in):: nk(klon) - real, intent(in):: tnk(klon), qnk(klon), gznk(klon) + real, intent(in):: tnk(:), qnk(:), gznk(:) ! (klon) real, intent(in):: t(klon, klev), qs(klon, klev), gz(klon, klev) real, intent(in):: p(klon, klev), h(klon, klev) - real, intent(in):: tv(klon, klev), lv(klon, klev) - real, intent(in):: pbase(klon), buoybase(klon), plcl(klon) + real, intent(in):: tv(klon, klev) + real, intent(in):: lv(:, :) ! (ncum, nl) + real, intent(in):: pbase(:), buoybase(:), plcl(:) ! (ncum) ! outputs: integer, intent(out):: inb(:) ! (ncum) @@ -39,10 +37,23 @@ real tp(klon, klev), tvp(klon, klev), clw(klon, klev) ! condensed water not removed from tvp - real hp(klon, klev), ep(klon, klev), sigp(klon, klev) + real hp(klon, klev), ep(klon, klev) real buoy(klon, klev) ! Local: + + integer ncum + + real, parameter:: pbcrit = 150. + ! critical cloud depth (mbar) beneath which the precipitation + ! efficiency is assumed to be zero + + real, parameter:: ptcrit = 500. + ! cloud depth (mbar) above which the precipitation efficiency is + ! assumed to be unity + + real, parameter:: dtovsh = - 0.2 ! dT for overshoot + integer i, k real tg, qg, ahg, alv, s, tc, es, denom real pden @@ -50,12 +61,13 @@ !--------------------------------------------------------------------- + ncum = size(icb) + ! SOME INITIALIZATIONS do k = 1, nl do i = 1, ncum - ep(i, k) = 0.0 - sigp(i, k) = spfac + ep(i, k) = 0. end do end do @@ -67,8 +79,8 @@ ! Calculate certain parcel quantities, including static energy do i = 1, ncum - ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) & - + qnk(i) * (lv0 - clmcpv * (tnk(i) - 273.15)) + gznk(i) + ah0(i) = (rcpd * (1. - qnk(i)) + rcw * qnk(i)) * tnk(i) & + + qnk(i) * (rlvtt - clmcpv * (tnk(i) - 273.15)) + gznk(i) end do ! Find lifted parcel quantities above cloud base @@ -78,20 +90,20 @@ if (k >= (icbs(i) + 1)) then tg = t(i, k) qg = qs(i, k) - alv = lv0 - clmcpv * (t(i, k) - 273.15) + alv = rlvtt - clmcpv * (t(i, k) - 273.15) ! First iteration. - s = cpd * (1. - qnk(i)) + cl * qnk(i) & - + alv * alv * qg / (rrv * t(i, k) * t(i, k)) + s = rcpd * (1. - qnk(i)) + rcw * qnk(i) & + + alv * alv * qg / (rv * t(i, k) * t(i, k)) s = 1. / s - ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) + ahg = rcpd * tg + (rcw - rcpd) * qnk(i) * tg + alv * qg + gz(i, k) tg = tg + s * (ah0(i) - ahg) tc = tg - 273.15 denom = 243.5 + tc - denom = MAX(denom, 1.0) + denom = MAX(denom, 1.) es = 6.112 * exp(17.67 * tc / denom) @@ -99,41 +111,39 @@ ! Second iteration. - ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) + ahg = rcpd * tg + (rcw - rcpd) * qnk(i) * tg + alv * qg + gz(i, k) tg = tg + s * (ah0(i) - ahg) tc = tg - 273.15 denom = 243.5 + tc - denom = MAX(denom, 1.0) + denom = MAX(denom, 1.) es = 6.112 * exp(17.67 * tc / denom) qg = eps * es / (p(i, k) - es * (1. - eps)) - alv = lv0 - clmcpv * (t(i, k) - 273.15) + alv = rlvtt - clmcpv * (t(i, k) - 273.15) ! no approximation: tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) & - / (cpd + (cl - cpd) * qnk(i)) + / (rcpd + (rcw - rcpd) * qnk(i)) clw(i, k) = qnk(i) - qg - clw(i, k) = max(0.0, clw(i, k)) + clw(i, k) = max(0., clw(i, k)) ! qg utilise au lieu du vrai mixing ratio rg: tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing endif end do end do - ! SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF - ! PRECIPITATION FALLING OUTSIDE OF CLOUD - ! THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) + ! SET THE PRECIPITATION EFFICIENCIES + ! It MAY BE a FUNCTION OF TP(I), P(I) AND CLW(I) do k = 1, nl do i = 1, ncum pden = ptcrit - pbcrit ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax - ep(i, k) = max(ep(i, k), 0.0) + ep(i, k) = max(ep(i, k), 0.) ep(i, k) = min(ep(i, k), epmax) - sigp(i, k) = spfac end do end do @@ -188,8 +198,8 @@ do k = minorig + 1, nl do i = 1, ncum - if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, nk(i)) & - + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k) + if (k >= icb(i) .and. k <= inb(i)) hp(i, k) = h(i, minorig) & + + (lv(i, k) + (rcpd - rcpv) * t(i, k)) * ep(i, k) * clw(i, k) end do end do