--- trunk/Sources/phylmd/CV30_routines/cv30_unsat.f 2016/03/21 15:36:26 186 +++ trunk/Sources/phylmd/CV30_routines/cv30_unsat.f 2016/04/14 15:15:56 190 @@ -4,77 +4,74 @@ contains - SUBROUTINE cv30_unsat(nloc, ncum, nd, na, icb, inb, t, rr, rs, gz, u, v, p, & - ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, mp, rp, & - up, vp, wt, water, evap, b) + SUBROUTINE cv30_unsat(icb, inb, t, q, qs, gz, u, v, p, ph, th, tv, lv, cpn, & + ep, sigp, clw, m, ment, elij, delt, plcl, mp, qp, up, vp, wt, water, & + evap, b) use cv30_param_m, only: nl, sigd - use cvflag, only: cvflag_grav - use cvthermo, only: cpd, ginv, grav + use cv_thermo_m, only: cpd, ginv, grav + USE dimphy, ONLY: klon, klev ! inputs: - integer, intent(in):: nloc, ncum, nd, na integer, intent(in):: icb(:), inb(:) ! (ncum) - real t(nloc, nd), rr(nloc, nd), rs(nloc, nd) - real gz(nloc, na) - real u(nloc, nd), v(nloc, nd) - real p(nloc, nd), ph(nloc, nd + 1) - real th(nloc, na) - real tv(nloc, na) - real lv(nloc, na) - real cpn(nloc, na) - real ep(nloc, na), sigp(nloc, na), clw(nloc, na) - real m(nloc, na), ment(nloc, na, na), elij(nloc, na, na) + real, intent(in):: t(:, :), q(:, :), qs(:, :) ! (klon, klev) + real, intent(in):: gz(:, :) ! (klon, klev) + real, intent(in):: u(:, :), v(:, :) ! (klon, klev) + real p(klon, klev), ph(klon, klev + 1) + real th(klon, klev) + real tv(klon, klev) + real lv(klon, klev) + real cpn(klon, klev) + real, intent(in):: ep(klon, klev), sigp(klon, klev), clw(klon, klev) + real m(klon, klev), ment(klon, klev, klev), elij(klon, klev, klev) real, intent(in):: delt - real plcl(nloc) + real plcl(klon) ! outputs: - real mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na) - real wt(nloc, na), water(nloc, na), evap(nloc, na) - real b(:, :) ! (nloc, na) + real, intent(out):: mp(klon, klev) + real, intent(out):: qp(:, :), up(:, :), vp(:, :) ! (ncum, nl) + real wt(klon, klev), water(klon, klev), evap(klon, klev) + real, intent(out):: b(:, :) ! (ncum, nl - 1) ! Local: + integer ncum integer i, j, il, num1 real tinv, delti real awat, afac, afac1, afac2, bfac real pr1, pr2, sigt, b6, c6, revap, tevap, delth real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf real ampmax - real lvcp(nloc, na) - real wdtrain(nloc) - logical lwork(nloc) + real lvcp(klon, klev) + real wdtrain(size(icb)) + logical lwork(size(icb)) !------------------------------------------------------ + ncum = size(icb) delti = 1. / delt tinv = 1. / 3. mp = 0. + b = 0. do i = 1, nl do il = 1, ncum - mp(il, i) = 0. - rp(il, i) = rr(il, i) + qp(il, i) = q(il, i) up(il, i) = u(il, i) vp(il, i) = v(il, i) wt(il, i) = 0.001 water(il, i) = 0. evap(il, i) = 0. - b(il, i) = 0. lvcp(il, i) = lv(il, i) / cpn(il, i) enddo enddo ! check whether ep(inb) = 0, if so, skip precipitating ! downdraft calculation + forall (il = 1:ncum) lwork(il) = ep(il, inb(il)) >= 1e-4 - do il = 1, ncum - lwork(il) = .TRUE. - if (ep(il, inb(il)) < 0.0001) lwork(il) = .FALSE. - enddo - - wdtrain(:ncum) = 0. + wdtrain = 0. - downdraft_loop: DO i = nl + 1, 1, - 1 + downdraft_loop: DO i = nl - 1, 1, - 1 num1 = 0 do il = 1, ncum @@ -89,11 +86,7 @@ do il = 1, ncum if (i <= inb(il) .and. lwork(il)) then - if (cvflag_grav) then - wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i) - else - wdtrain(il) = 10. * ep(il, i) * m(il, i) * clw(il, i) - endif + wdtrain(il) = grav * ep(il, i) * m(il, i) * clw(il, i) endif enddo @@ -102,51 +95,48 @@ do il = 1, ncum if (i <= inb(il) .and. lwork(il)) then awat = elij(il, j, i) - (1. - ep(il, i)) * clw(il, i) - awat = amax1(awat, 0.) - if (cvflag_grav) then - wdtrain(il) = wdtrain(il) + grav * awat & - * ment(il, j, i) - else - wdtrain(il) = wdtrain(il) + 10. * awat * ment(il, j, i) - endif + awat = max(awat, 0.) + wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i) endif enddo end do endif ! find rain water and evaporation using provisional - ! estimates of rp(i)and rp(i - 1) + ! estimates of qp(i) and qp(i - 1) do il = 1, ncum if (i <= inb(il) .and. lwork(il)) then wt(il, i) = 45. if (i < inb(il)) then - rp(il, i) = rp(il, i + 1) + (cpd * (t(il, i + 1) & + qp(il, i) = qp(il, i + 1) + (cpd * (t(il, i + 1) & - t(il, i)) + gz(il, i + 1) - gz(il, i)) / lv(il, i) - rp(il, i) = 0.5 * (rp(il, i) + rr(il, i)) + qp(il, i) = 0.5 * (qp(il, i) + q(il, i)) endif - rp(il, i) = amax1(rp(il, i), 0.) - rp(il, i) = amin1(rp(il, i), rs(il, i)) - rp(il, inb(il)) = rr(il, inb(il)) + + qp(il, i) = max(qp(il, i), 0.) + qp(il, i) = min(qp(il, i), qs(il, i)) + qp(il, inb(il)) = q(il, inb(il)) if (i == 1) then - afac = p(il, 1) * (rs(il, 1) - rp(il, 1)) & - / (1e4 + 2000. * p(il, 1) * rs(il, 1)) + afac = p(il, 1) * (qs(il, 1) - qp(il, 1)) & + / (1e4 + 2000. * p(il, 1) * qs(il, 1)) else - rp(il, i - 1) = rp(il, i) + (cpd * (t(il, i) & + qp(il, i - 1) = qp(il, i) + (cpd * (t(il, i) & - t(il, i - 1)) + gz(il, i) - gz(il, i - 1)) / lv(il, i) - rp(il, i - 1) = 0.5 * (rp(il, i - 1) + rr(il, i - 1)) - rp(il, i - 1) = amin1(rp(il, i - 1), rs(il, i - 1)) - rp(il, i - 1) = amax1(rp(il, i - 1), 0.) - afac1 = p(il, i) * (rs(il, i) - rp(il, i)) & - / (1e4 + 2000. * p(il, i) * rs(il, i)) - afac2 = p(il, i - 1) * (rs(il, i - 1) - rp(il, i - 1)) & - / (1e4 + 2000. * p(il, i - 1) * rs(il, i - 1)) + qp(il, i - 1) = 0.5 * (qp(il, i - 1) + q(il, i - 1)) + qp(il, i - 1) = min(qp(il, i - 1), qs(il, i - 1)) + qp(il, i - 1) = max(qp(il, i - 1), 0.) + afac1 = p(il, i) * (qs(il, i) - qp(il, i)) & + / (1e4 + 2000. * p(il, i) * qs(il, i)) + afac2 = p(il, i - 1) * (qs(il, i - 1) - qp(il, i - 1)) & + / (1e4 + 2000. * p(il, i - 1) * qs(il, i - 1)) afac = 0.5 * (afac1 + afac2) endif - if (i == inb(il))afac = 0. - afac = amax1(afac, 0.) + + if (i == inb(il)) afac = 0. + afac = max(afac, 0.) bfac = 1. / (sigd * wt(il, i)) ! prise en compte de la variation progressive de sigt dans @@ -180,26 +170,22 @@ ! hydrostatic approximation if (i /= 1) then - tevap = amax1(0., evap(il, i)) - delth = amax1(0.001, (th(il, i) - th(il, i - 1))) - if (cvflag_grav) then - mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap & - * (p(il, i - 1) - p(il, i)) / delth - else - mp(il, i) = 10. * lvcp(il, i) * sigd * tevap & - * (p(il, i - 1) - p(il, i)) / delth - endif - - ! if hydrostatic assumption fails, - ! solve cubic difference equation for downdraft theta - ! and mass flux from two simultaneous differential eqns + tevap = max(0., evap(il, i)) + delth = max(0.001, (th(il, i) - th(il, i - 1))) + mp(il, i) = 100. * ginv * lvcp(il, i) * sigd * tevap & + * (p(il, i - 1) - p(il, i)) / delth + + ! If hydrostatic assumption fails, solve cubic + ! difference equation for downdraft theta and mass + ! flux from two simultaneous differential equations amfac = sigd * sigd * 70. * ph(il, i) & * (p(il, i - 1) - p(il, i)) & * (th(il, i) - th(il, i - 1)) / (tv(il, i) * th(il, i)) amp2 = abs(mp(il, i + 1) * mp(il, i + 1) - mp(il, i) & * mp(il, i)) - if (amp2 > (0.1 * amfac)) then + + if (amp2 > 0.1 * amfac) then xf = 100. * sigd * sigd * sigd * (ph(il, i) & - ph(il, i + 1)) tf = b(il, i) - 5. * (th(il, i) - th(il, i - 1)) & @@ -209,13 +195,14 @@ * mp(il, i + 1) * xf * tf + 50. * (p(il, i - 1) & - p(il, i)) * xf * tevap fac2 = 1. - if (bf < 0.)fac2 = - 1. + if (bf < 0.) fac2 = - 1. bf = abs(bf) ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv + if (ur >= 0.) then sru = sqrt(ur) fac = 1. - if ((0.5 * bf - sru) < 0.)fac = - 1. + if ((0.5 * bf - sru) < 0.) fac = - 1. mp(il, i) = mp(il, i + 1) * tinv & + (0.5 * bf + sru)**tinv & + fac * (abs(0.5 * bf - sru))**tinv @@ -225,61 +212,44 @@ mp(il, i) = mp(il, i + 1) * tinv + 2. & * sqrt(af * tinv) * cos(d * tinv) endif - mp(il, i) = amax1(0., mp(il, i)) - if (cvflag_grav) then - ! Il y a vraisemblablement une erreur dans la - ! ligne 2 suivante: il faut diviser par (mp(il, - ! i) * sigd * grav) et non par (mp(il, i) + sigd - ! * 0.1). Et il faut bien revoir les facteurs - ! 100. - b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) & - - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) & - - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) & - / (lvcp(il, i) * sigd * th(il, i)) - else - b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) & - - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) & - - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) & - / (lvcp(il, i) * sigd * th(il, i)) - endif - b(il, i - 1) = amax1(b(il, i - 1), 0.) + mp(il, i) = max(0., mp(il, i)) + + ! Il y a vraisemblablement une erreur dans la + ! ligne suivante : il faut diviser par (mp(il, + ! i) * sigd * grav) et non par (mp(il, i) + sigd + ! * 0.1). Et il faut bien revoir les facteurs + ! 100. + b(il, i - 1) = b(il, i) + 100. * (p(il, i - 1) & + - p(il, i)) * tevap / (mp(il, i) + sigd * 0.1) & + - 10. * (th(il, i) - th(il, i - 1)) * t(il, i) & + / (lvcp(il, i) * sigd * th(il, i)) + b(il, i - 1) = max(b(il, i - 1), 0.) endif ! limit magnitude of mp(i) to meet cfl condition - ampmax = 2. * (ph(il, i) - ph(il, i + 1)) * delti amp2 = 2. * (ph(il, i - 1) - ph(il, i)) * delti - ampmax = amin1(ampmax, amp2) - mp(il, i) = amin1(mp(il, i), ampmax) + ampmax = min(ampmax, amp2) + mp(il, i) = min(mp(il, i), ampmax) ! force mp to decrease linearly to zero ! between cloud base and the surface - - if (p(il, i) > p(il, icb(il))) then - mp(il, i) = mp(il, icb(il)) * (p(il, 1) - p(il, i)) & - / (p(il, 1) - p(il, icb(il))) - endif + if (p(il, i) > p(il, icb(il))) mp(il, i) = mp(il, icb(il)) & + * (p(il, 1) - p(il, i)) / (p(il, 1) - p(il, icb(il))) endif ! i == 1 ! find mixing ratio of precipitating downdraft if (i /= inb(il)) then - rp(il, i) = rr(il, i) + qp(il, i) = q(il, i) if (mp(il, i) > mp(il, i + 1)) then - if (cvflag_grav) then - rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) & - * (mp(il, i) - mp(il, i + 1)) + 100. * ginv & - * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) & - * (evap(il, i + 1) + evap(il, i)) - else - rp(il, i) = rp(il, i + 1) * mp(il, i + 1) + rr(il, i) & - * (mp(il, i) - mp(il, i + 1)) + 5. * sigd & - * (ph(il, i) - ph(il, i + 1)) & - * (evap(il, i + 1) + evap(il, i)) - endif - rp(il, i) = rp(il, i) / mp(il, i) + qp(il, i) = qp(il, i + 1) * mp(il, i + 1) + q(il, i) & + * (mp(il, i) - mp(il, i + 1)) + 100. * ginv & + * 0.5 * sigd * (ph(il, i) - ph(il, i + 1)) & + * (evap(il, i + 1) + evap(il, i)) + qp(il, i) = qp(il, i) / mp(il, i) up(il, i) = up(il, i + 1) * mp(il, i + 1) + u(il, i) & * (mp(il, i) - mp(il, i + 1)) up(il, i) = up(il, i) / mp(il, i) @@ -288,24 +258,16 @@ vp(il, i) = vp(il, i) / mp(il, i) else if (mp(il, i + 1) > 1e-16) then - if (cvflag_grav) then - rp(il, i) = rp(il, i + 1) & - + 100. * ginv * 0.5 * sigd * (ph(il, i) & - - ph(il, i + 1)) & - * (evap(il, i + 1) + evap(il, i)) & - / mp(il, i + 1) - else - rp(il, i) = rp(il, i + 1) & - + 5. * sigd * (ph(il, i) - ph(il, i + 1)) & - * (evap(il, i + 1) + evap(il, i)) & - / mp(il, i + 1) - endif + qp(il, i) = qp(il, i + 1) + 100. * ginv * 0.5 * sigd & + * (ph(il, i) - ph(il, i + 1)) & + * (evap(il, i + 1) + evap(il, i)) / mp(il, i + 1) up(il, i) = up(il, i + 1) vp(il, i) = vp(il, i + 1) endif endif - rp(il, i) = amin1(rp(il, i), rs(il, i)) - rp(il, i) = amax1(rp(il, i), 0.) + + qp(il, i) = min(qp(il, i), qs(il, i)) + qp(il, i) = max(qp(il, i), 0.) endif endif end do