--- trunk/phylmd/CV3_routines/cv3_yield.f 2014/04/25 14:58:31 97 +++ trunk/phylmd/CV30_routines/cv30_yield.f 2018/02/05 10:39:38 254 @@ -1,708 +1,583 @@ -module cv3_yield_m +module cv30_yield_m implicit none contains - SUBROUTINE cv3_yield(nloc,ncum,nd,na & - ,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,qent,uent,vent,nent,elij,sig & - ,tv,tvp & - ,iflag,precip,VPrecip,ft,fr,fu,fv & - ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd) - use conema3_m - use cv3_param_m - use cvthermo - use cvflag + SUBROUTINE cv30_yield(icb, inb, delt, t, rr, u, v, gz, p, ph, h, hp, lv, & + 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, ma, mike, tls, tps, qcondc) + + ! Tendencies, precipitation, variables of interface with other + ! processes, etc. + + use conf_phys_m, only: iflag_clw + 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):: ncum,nd,na,nloc - integer icb(nloc), inb(nloc) - real, intent(in):: delt - real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd) - real sig(nloc,nd) - real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na) - real th(nloc,na), p(nloc,nd), tp(nloc,na) - real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na) - real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na) - real vp(nloc,na), wt(nloc,nd) - real water(nloc,na), evap(nloc,na), b(nloc,na) - real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na) - !ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na) - real vent(nloc,na,na), elij(nloc,na,na) - integer nent(nloc,na) - real tv(nloc,nd), tvp(nloc,nd) - ! input/output: - integer iflag(nloc) + 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, 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, 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, intent(in):: mp(:, :) ! (ncum, nl) Mass flux of the + ! unsaturated downdraft, defined positive downward, in kg m-2 + ! s-1. M_p in Emanuel (1991 928). + + 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) + real vent(klon, klev, klev) + integer nent(klon, klev) + real elij(klon, klev, klev) + real sig(klon, klev) + real tv(klon, klev), tvp(klon, klev) ! outputs: - real precip(nloc) - real VPrecip(nloc,nd+1) - real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd) - real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd) - real dnwd0(nloc,nd), mike(nloc,nd) - real tls(nloc,nd), tps(nloc,nd) - real qcondc(nloc,nd) ! cld - real wd(nloc) ! gust - - ! local variables: - integer i,k,il,n,j,num1 - real rat, awat, delti - real ax, bx, cx, dx, ex + 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) + real upwd(klon, klev), dnwd(klon, klev) + real ma(klon, klev) + real mike(klon, klev) + real tls(klon, klev), tps(klon, klev) + real qcondc(klon, klev) + + ! Local: + real, parameter:: delta = 0.01 ! interface cloud parameterization + integer ncum + integer i, k, il, n, j + real awat, delti + real ax, bx, cx, dx real cpinv, rdcp, dpinv - real lvcp(nloc,na), mke(nloc,na) - real am(nloc), work(nloc), ad(nloc), amp1(nloc) -!!! real up1(nloc), dn1(nloc) - real up1(nloc,nd,nd), dn1(nloc,nd,nd) - real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc) - real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld - real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld - + real lvcp(klon, klev) + real am(klon), work(klon), ad(klon), amp1(klon) + real up1(klon, klev, klev), dn1(klon, klev, klev) + real asum(klon), bsum(klon), csum(klon), dsum(klon) + real qcond(klon, klev), nqcond(klon, klev), wa(klon, klev) + real siga(klon, klev), sax(klon, klev), mac(klon, klev) !------------------------------------------------------------- + ncum = size(icb) + iflag = 0 + ! initialization: - delti = 1.0/delt + delti = 1.0 / delt - do il=1,ncum - precip(il)=0.0 - wd(il)=0.0 ! gust - VPrecip(il,nd+1)=0. + do il = 1, ncum + precip(il) = 0.0 + VPrecip(il, klev + 1) = 0. enddo - do i=1,nd - do il=1,ncum - VPrecip(il,i)=0.0 - ft(il,i)=0.0 - fr(il,i)=0.0 - fu(il,i)=0.0 - fv(il,i)=0.0 - qcondc(il,i)=0.0 ! cld - qcond(il,i)=0.0 ! cld - nqcond(il,i)=0.0 ! cld + do i = 1, klev + do il = 1, ncum + VPrecip(il, i) = 0.0 + ft(il, i) = 0.0 + fr(il, i) = 0.0 + fu(il, i) = 0.0 + fv(il, i) = 0.0 + qcondc(il, i) = 0.0 + qcond(il, i) = 0.0 + nqcond(il, i) = 0.0 enddo enddo - - do i=1,nl - do il=1,ncum - lvcp(il,i)=lv(il,i)/cpn(il,i) + do i = 1, nl + do il = 1, ncum + lvcp(il, i) = lv(il, i) / cpn(il, i) enddo enddo + ! calculate surface precipitation in mm / day - ! - ! *** calculate surface precipitation in mm/day *** - ! - do il=1,ncum - if(ep(il,inb(il)).ge.0.0001)then - if (cvflag_grav) then - precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav) - else - precip(il)=wt(il,1)*sigd*water(il,1)*8640. - endif - endif + do il = 1, ncum + if (ep(il, inb(il)) >= 1e-4) precip(il) = wt(il, 1) * sigd & + * water(il, 1) * 86400. * 1000. / (rowl * rg) enddo - ! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s === - ! + ! CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg / m2 / s === + ! MAF rajout pour lessivage - do k=1,nl - do il=1,ncum - if (k.le.inb(il)) then - if (cvflag_grav) then - VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav - else - VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10. - endif - endif + do k = 1, nl - 1 + do il = 1, ncum + if (k <= inb(il)) VPrecip(il, k) = wt(il, k) * sigd * water(il, k) & + / rg end do end do - ! - ! - ! - ! *** calculate tendencies of lowest level potential temperature *** - ! *** and mixing ratio *** - ! - do il=1,ncum - work(il)=1.0/(ph(il,1)-ph(il,2)) - am(il)=0.0 + + ! calculate tendencies of lowest level potential temperature + ! and mixing ratio + + do il = 1, ncum + work(il) = 1.0 / (ph(il, 1) - ph(il, 2)) + am(il) = 0.0 + enddo + + do k = 2, nl + do il = 1, ncum + if (k <= inb(il)) am(il) = am(il) + m(il, k) + enddo enddo - do k=2,nl - do il=1,ncum - if (k.le.inb(il)) then - am(il)=am(il)+m(il,k) + do il = 1, ncum + if (0.01 * rg * work(il) * am(il) >= delti) iflag(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) + + ! jyg Correction pour mieux conserver l'eau (conformite avec CONVECT4.3) + ! (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 * rg * am(il) * (rr(il, 2) - rr(il, 1)) & + * work(il) + + 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 * rg * work(il) * (mp(il, 2) & + * (vp(il, 2) - v(il, 1)) + am(il) * (v(il, 2) - v(il, 1))) + enddo + + do j = 2, nl + do il = 1, ncum + if (j <= inb(il)) then + 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 * rg * work(il) * ment(il, j, 1) & + * (uent(il, j, 1) - u(il, 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 - do il=1,ncum + ! calculate tendencies of potential temperature and mixing ratio + ! at levels above the lowest level - ! convect3 if((0.1*dpinv*am).ge.delti)iflag(il)=4 - if (cvflag_grav) then - if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect - 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)) - else - if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect - ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1) & - +(gz(il,2)-gz(il,1))/cpn(il,1)) - endif - - ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2)) - - if (cvflag_grav) then - ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2) & - *t(il,1)*b(il,1)*work(il) - else - ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il) - endif - - 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) - - if (cvflag_grav) then - !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))*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))*work(il) - - fu(il,1)=fu(il,1)+0.01*grav*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)*(vp(il,2)-v(il,1)) & - +am(il)*(v(il,2)-v(il,1))) - else ! cvflag_grav - fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) & - +sigd*0.5*(evap(il,1)+evap(il,2)) - fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il) - fu(il,1)=fu(il,1)+0.1*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.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) & - +am(il)*(v(il,2)-v(il,1))) - endif ! cvflag_grav - - enddo ! il - - do j=2,nl - do il=1,ncum - if (j.le.inb(il)) then - if (cvflag_grav) then - fr(il,1)=fr(il,1) & - +0.01*grav*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)*(uent(il,j,1)-u(il,1)) - fv(il,1)=fv(il,1) & - +0.01*grav*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1)) - else ! cvflag_grav - fr(il,1)=fr(il,1) & - +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1)) - fu(il,1)=fu(il,1) & - +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1)) - fv(il,1)=fv(il,1) & - +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1)) - endif ! cvflag_grav - endif ! j - enddo - enddo - - ! - ! *** 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 *** - ! - - do i=2,nl+1 ! newvecto: mettre nl au lieu nl+1? - - num1=0 - do il=1,ncum - if(i.le.inb(il))num1=num1+1 - enddo - if(num1.le.0) cycle - - call zilch(amp1,ncum) - call zilch(ad,ncum) - - do k=i+1,nl+1 - do il=1,ncum - if (i.le.inb(il) .and. k.le.(inb(il)+1)) then - amp1(il)=amp1(il)+m(il,k) - endif - end do - end do + ! first find the net saturated updraft and downdraft mass fluxes + ! through each level - do k=1,i - do j=i+1,nl+1 - do il=1,ncum - if (i.le.inb(il) .and. j.le.(inb(il)+1)) then - amp1(il)=amp1(il)+ment(il,k,j) + loop_i: do i = 2, nl - 1 + if (any(inb >= i)) then + amp1(:ncum) = 0. + ad(:ncum) = 0. + + do k = i + 1, nl + 1 + do il = 1, ncum + if (i <= inb(il) .and. k <= (inb(il) + 1)) then + amp1(il) = amp1(il) + m(il, k) endif end do end do - end do - do k=1,i-1 - do j=i,nl+1 ! newvecto: nl au lieu nl+1? - do il=1,ncum - if (i.le.inb(il) .and. j.le.inb(il)) then - ad(il)=ad(il)+ment(il,j,k) - endif + do k = 1, i + do j = i + 1, nl + 1 + do il = 1, ncum + if (i <= inb(il) .and. j <= (inb(il) + 1)) then + amp1(il) = amp1(il) + ment(il, k, j) + endif + end do end do end do - end do - do il=1,ncum - if (i.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) - - ! convect3 if((0.1*dpinv*amp1).ge.delti)iflag(il)=4 - if (cvflag_grav) then - if((0.01*grav*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto - else - if((0.1*dpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto - endif - - if (cvflag_grav) then - ft(il,i)=0.01*grav*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 - else ! cvflag_grav - ft(il,i)=0.1*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.09*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.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i) & - +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv - endif ! cvflag_grav - - - 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 - - if (cvflag_grav) then - fr(il,i)=0.01*grav*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)*(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)*(v(il,i+1)-v(il,i)) & - -ad(il)*(v(il,i)-v(il,i-1))) - else ! cvflag_grav - fr(il,i)=0.1*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.1*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.1*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) & - -ad(il)*(v(il,i)-v(il,i-1))) - endif ! cvflag_grav + do k = 1, i - 1 + do j = i, nl + 1 ! newvecto: nl au lieu nl + 1? + do il = 1, ncum + if (i <= inb(il) .and. j <= inb(il)) then + ad(il) = ad(il) + ment(il, j, k) + endif + end do + end do + end do - endif ! i - end do + do il = 1, ncum + if (i <= inb(il)) then + dpinv = 1.0 / (ph(il, i) - ph(il, i + 1)) + cpinv = 1.0 / cpn(il, i) + + if (0.01 * rg * dpinv * amp1(il) >= delti) iflag(il) = 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)) - 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 * 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 * rg * dpinv * (amp1(il) & + * (v(il, i + 1) - v(il, i)) - ad(il) * (v(il, i) & + - v(il, i - 1))) + endif + end do - do k=1,i-1 - do il=1,ncum - if (i.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) - - awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i) - awat=amax1(awat,0.0) - - if (cvflag_grav) then - fr(il,i)=fr(il,i) & - +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i)) - fu(il,i)=fu(il,i) & - +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) - fv(il,i)=fv(il,i) & - +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) - else ! cvflag_grav - fr(il,i)=fr(il,i) & - +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i)) - fu(il,i)=fu(il,i) & - +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) - fv(il,i)=fv(il,i) & - +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) - endif ! cvflag_grav - - ! (saturated updrafts resulting from mixing) ! cld - qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld - nqcond(il,i)=nqcond(il,i)+1. ! cld - endif ! i + do k = 1, i - 1 + do il = 1, ncum + if (i <= inb(il)) then + dpinv = 1.0 / (ph(il, i) - ph(il, i + 1)) + cpinv = 1.0 / cpn(il, i) + + 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 * rg * dpinv & + * ment(il, k, i) * (qent(il, k, i) - awat - rr(il, i)) + 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 * rg * dpinv & + * ment(il, k, i) * (vent(il, k, i) - v(il, i)) + + ! (saturated updrafts resulting from mixing) + qcond(il, i) = qcond(il, i) + (elij(il, k, i) - awat) + nqcond(il, i) = nqcond(il, i) + 1. + endif ! i + end do end do - end do - do k=i,nl+1 - do il=1,ncum - if (i.le.inb(il) .and. k.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) - - if (cvflag_grav) then - fr(il,i)=fr(il,i) & - +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i)) - fu(il,i)=fu(il,i) & - +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) - fv(il,i)=fv(il,i) & - +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) - else ! cvflag_grav - fr(il,i)=fr(il,i) & - +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i)) - fu(il,i)=fu(il,i) & - +0.1*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i)) - fv(il,i)=fv(il,i) & - +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i)) - endif ! cvflag_grav - endif ! i and k + do k = i, nl + 1 + do il = 1, ncum + if (i <= inb(il) .and. k <= inb(il)) then + dpinv = 1.0 / (ph(il, i) - ph(il, i + 1)) + cpinv = 1.0 / cpn(il, i) + + 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 * rg * dpinv & + * ment(il, k, i) * (uent(il, k, i) - u(il, i)) + fv(il, i) = fv(il, i) + 0.01 * rg * dpinv & + * ment(il, k, i) * (vent(il, k, i) - v(il, i)) + endif + end do end do - end do - do il=1,ncum - if (i.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) + do il = 1, ncum + if (i <= inb(il)) then + dpinv = 1.0 / (ph(il, i) - ph(il, i + 1)) + cpinv = 1.0 / cpn(il, i) - if (cvflag_grav) then ! 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)-rr(il,i-1)))*dpinv - - fu(il,i)=fu(il,i)+0.01*grav*(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)*(vp(il,i+1)-v(il,i)) & - -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv - else ! cvflag_grav - fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1)) & - +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) & - *(rp(il,i)-rr(il,i-1)))*dpinv - fu(il,i)=fu(il,i)+0.1*(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.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) & - -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv - endif ! cvflag_grav - - endif ! i - end do + fr(il, i) = fr(il, i) + 0.5 * sigd * (evap(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 * 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 * rg * (mp(il, i + 1) & + * (vp(il, i + 1) - v(il, i)) - mp(il, i) * (vp(il, i) & + - v(il, i - 1))) * dpinv + endif + end do - ! sb: interface with the cloud parameterization: ! cld + ! sb: interface with the cloud parameterization: - do k=i+1,nl - do il=1,ncum - if (k.le.inb(il) .and. i.le.inb(il)) then ! cld - ! (saturated downdrafts resulting from mixing) ! cld - qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld - nqcond(il,i)=nqcond(il,i)+1. ! cld - endif ! cld - enddo ! cld - enddo ! cld - - ! (particular case: no detraining level is found) ! cld - do il=1,ncum ! cld - if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld - qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld - nqcond(il,i)=nqcond(il,i)+1. ! cld - endif ! cld - enddo ! cld - - do il=1,ncum ! cld - if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld - qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld - endif ! cld - enddo + do k = i + 1, nl + do il = 1, ncum + if (k <= inb(il) .and. i <= inb(il)) then + ! (saturated downdrafts resulting from mixing) + qcond(il, i) = qcond(il, i) + elij(il, k, i) + nqcond(il, i) = nqcond(il, i) + 1. + endif + enddo + enddo - end do + ! (particular case: no detraining level is found) + do il = 1, ncum + if (i <= inb(il) .and. nent(il, i) == 0) then + qcond(il, i) = qcond(il, i) + (1. - ep(il, i)) * clw(il, i) + nqcond(il, i) = nqcond(il, i) + 1. + endif + enddo + do il = 1, ncum + if (i <= inb(il) .and. nqcond(il, i) /= 0.) then + qcond(il, i) = qcond(il, i) / nqcond(il, i) + endif + enddo + 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 *** - ! - 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) & - *(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 - ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il)) & - *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1) & - *(ph(il,inb(il)-1)-ph(il,inb(il)))) - - bx=0.1*ment(il,inb(il),inb(il))*(qent(il,inb(il),inb(il)) & - -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) - fr(il,inb(il))=fr(il,inb(il))-bx - fr(il,inb(il)-1)=fr(il,inb(il)-1) & - +bx*(ph(il,inb(il))-ph(il,inb(il)+1)) & - /(ph(il,inb(il)-1)-ph(il,inb(il))) - - cx=0.1*ment(il,inb(il),inb(il))*(uent(il,inb(il),inb(il)) & - -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) - fu(il,inb(il))=fu(il,inb(il))-cx - fu(il,inb(il)-1)=fu(il,inb(il)-1) & - +cx*(ph(il,inb(il))-ph(il,inb(il)+1)) & - /(ph(il,inb(il)-1)-ph(il,inb(il))) - - dx=0.1*ment(il,inb(il),inb(il))*(vent(il,inb(il),inb(il)) & - -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) - fv(il,inb(il))=fv(il,inb(il))-dx - fv(il,inb(il)-1)=fv(il,inb(il)-1) & - +dx*(ph(il,inb(il))-ph(il,inb(il)+1)) & - /(ph(il,inb(il)-1)-ph(il,inb(il))) + ! 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)) * (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 + ft(il, inb(il) - 1) = ft(il, inb(il) - 1) + ax * cpn(il, inb(il)) & + * (ph(il, inb(il)) - ph(il, inb(il) + 1)) / (cpn(il, inb(il) - 1) & + * (ph(il, inb(il) - 1) - ph(il, inb(il)))) + + bx = 0.1 * ment(il, inb(il), inb(il)) * (qent(il, inb(il), inb(il)) & + - rr(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1)) + fr(il, inb(il)) = fr(il, inb(il)) - bx + fr(il, inb(il) - 1) = fr(il, inb(il) - 1) & + + bx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) & + / (ph(il, inb(il) - 1) - ph(il, inb(il))) + + cx = 0.1 * ment(il, inb(il), inb(il)) * (uent(il, inb(il), inb(il)) & + - u(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1)) + fu(il, inb(il)) = fu(il, inb(il)) - cx + fu(il, inb(il) - 1) = fu(il, inb(il) - 1) & + + cx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) & + / (ph(il, inb(il) - 1) - ph(il, inb(il))) + + dx = 0.1 * ment(il, inb(il), inb(il)) * (vent(il, inb(il), inb(il)) & + - v(il, inb(il))) / (ph(il, inb(il)) - ph(il, inb(il) + 1)) + fv(il, inb(il)) = fv(il, inb(il)) - dx + fv(il, inb(il) - 1) = fv(il, inb(il) - 1) & + + dx * (ph(il, inb(il)) - ph(il, inb(il) + 1)) & + / (ph(il, inb(il) - 1) - ph(il, inb(il))) end do - ! - ! *** homoginize tendencies below cloud base *** - ! - ! - do il=1,ncum - asum(il)=0.0 - bsum(il)=0.0 - csum(il)=0.0 - dsum(il)=0.0 - enddo - - do i=1,nl - do il=1,ncum - if (i.le.(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)*(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)-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) + ! homoginize tendencies below cloud base + + do il = 1, ncum + asum(il) = 0.0 + bsum(il) = 0.0 + csum(il) = 0.0 + dsum(il) = 0.0 + enddo + + do i = 1, nl + 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) + (rcw - rcpd) & + * (t(il, i) - t(il, 1))) * (ph(il, i) - ph(il, i + 1)) + 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) endif enddo enddo -!!!! do 700 i=1,icb(il)-1 - do i=1,nl - do il=1,ncum - if (i.le.(icb(il)-1)) then - ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il)) - fr(il,i)=bsum(il)/csum(il) + do i = 1, nl + do il = 1, ncum + if (i <= (icb(il) - 1)) then + ft(il, i) = asum(il) * t(il, i) / (th(il, i) * dsum(il)) + fr(il, i) = bsum(il) / csum(il) endif enddo enddo - ! - ! *** reset counter and return *** - ! - do il=1,ncum - sig(il,nd)=2.0 - enddo - + ! reset counter and return - do i=1,nd - do il=1,ncum - upwd(il,i)=0.0 - dnwd(il,i)=0.0 - enddo + do il = 1, ncum + sig(il, klev) = 2.0 enddo - do i=1,nl - do il=1,ncum - dnwd0(il,i)=-mp(il,i) + do i = 1, klev + do il = 1, ncum + upwd(il, i) = 0.0 + dnwd(il, i) = 0.0 enddo enddo - do i=nl+1,nd - do il=1,ncum - dnwd0(il,i)=0. - enddo - enddo - - do i=1,nl - do il=1,ncum - if (i.ge.icb(il) .and. i.le.inb(il)) then - upwd(il,i)=0.0 - dnwd(il,i)=0.0 + do i = 1, nl + do il = 1, ncum + if (i >= icb(il) .and. i <= inb(il)) then + upwd(il, i) = 0.0 + dnwd(il, i) = 0.0 endif enddo enddo - do i=1,nl - do k=1,nl - do il=1,ncum - up1(il,k,i)=0.0 - dn1(il,k,i)=0.0 + do i = 1, nl + do k = 1, nl + do il = 1, ncum + up1(il, k, i) = 0.0 + dn1(il, k, i) = 0.0 enddo enddo enddo - do i=1,nl - do k=i,nl - do n=1,i-1 - do il=1,ncum - if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then - up1(il,k,i)=up1(il,k,i)+ment(il,n,k) - dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n) + do i = 1, nl + do k = i, nl + do n = 1, i - 1 + do il = 1, ncum + if (i >= icb(il).and.i <= inb(il).and.k <= inb(il)) then + up1(il, k, i) = up1(il, k, i) + ment(il, n, k) + dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n) endif enddo enddo enddo enddo - do i=2,nl - do k=i,nl - do il=1,ncum - !test if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then - if (i.le.inb(il).and.k.le.inb(il)) then - upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i) - dnwd(il,i)=dnwd(il,i)+dn1(il,k,i) + do i = 2, nl + do k = i, nl + do il = 1, ncum + if (i <= inb(il).and.k <= inb(il)) then + upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i) + dnwd(il, i) = dnwd(il, i) + dn1(il, k, i) endif enddo enddo enddo + ! D\'etermination de la variation de flux ascendant entre + ! deux niveaux non dilu\'es mike - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! determination de la variation de flux ascendant entre - ! deux niveau non dilue mike - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - do i=1,nl - do il=1,ncum - mike(il,i)=m(il,i) + do i = 1, nl + do il = 1, ncum + mike(il, i) = m(il, i) enddo enddo - do i=nl+1,nd - do il=1,ncum - mike(il,i)=0. + do i = nl + 1, klev + do il = 1, ncum + mike(il, i) = 0. enddo enddo - do i=1,nd - do il=1,ncum - ma(il,i)=0 + do i = 1, klev + do il = 1, ncum + ma(il, i) = 0 enddo enddo - do i=1,nl - do j=i,nl - do il=1,ncum - ma(il,i)=ma(il,i)+m(il,j) + do i = 1, nl + do j = i, nl + do il = 1, ncum + ma(il, i) = ma(il, i) + m(il, j) enddo enddo enddo - do i=nl+1,nd - do il=1,ncum - ma(il,i)=0. + do i = nl + 1, klev + do il = 1, ncum + ma(il, i) = 0. enddo enddo - do i=1,nl - do il=1,ncum - if (i.le.(icb(il)-1)) then - ma(il,i)=0 + do i = 1, nl + do il = 1, ncum + if (i <= (icb(il) - 1)) then + ma(il, i) = 0 endif enddo enddo - !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ! icb represente de niveau ou se trouve la - ! base du nuage , et inb le top du nuage - !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - do i=1,nd - do il=1,ncum - mke(il,i)=upwd(il,i)+dnwd(il,i) + ! 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 = (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 + + do i = 1, klev + do il = 1, ncum + mac(il, i) = 0.0 + wa(il, i) = 0.0 + siga(il, i) = 0.0 + sax(il, i) = 0.0 + enddo + enddo + + do i = minorig, nl + do k = i + 1, nl + 1 + do il = 1, ncum + if (i <= inb(il) .and. k <= (inb(il) + 1)) then + mac(il, i) = mac(il, i) + m(il, k) + endif + enddo enddo enddo - do i=1,nd - DO il=1,ncum - rdcp=(rrd*(1.-rr(il,i))-rr(il,i)*rrv) & - /(cpd*(1.-rr(il,i))+rr(il,i)*cpv) - tls(il,i)=t(il,i)*(1000.0/p(il,i))**rdcp - tps(il,i)=tp(il,i) - end DO + do i = 1, nl + do j = 1, i + do il = 1, ncum + if (i >= icb(il) .and. i <= (inb(il) - 1) & + .and. j >= icb(il)) then + 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 + enddo enddo - ! - ! *** diagnose the in-cloud mixing ratio *** ! cld - ! *** of condensed water *** ! cld - ! ! cld - - do i=1,nd ! cld - do il=1,ncum ! cld - mac(il,i)=0.0 ! cld - wa(il,i)=0.0 ! cld - siga(il,i)=0.0 ! cld - sax(il,i)=0.0 ! cld - enddo ! cld - enddo ! cld - - do i=minorig, nl ! cld - do k=i+1,nl+1 ! cld - do il=1,ncum ! cld - if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld - mac(il,i)=mac(il,i)+m(il,k) ! cld - endif ! cld - enddo ! cld - enddo ! cld - enddo ! cld - - do i=1,nl ! cld - do j=1,i ! cld - do il=1,ncum ! cld - if (i.ge.icb(il) .and. i.le.(inb(il)-1) & - .and. j.ge.icb(il) ) then ! cld - sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) & - *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld - endif ! cld - enddo ! cld - enddo ! cld - enddo ! cld - - do i=1,nl ! cld - do il=1,ncum ! cld - if (i.ge.icb(il) .and. i.le.(inb(il)-1) & - .and. sax(il,i).gt.0.0 ) then ! cld - wa(il,i)=sqrt(2.*sax(il,i)) ! cld - endif ! cld - enddo ! cld - enddo ! cld - - do i=1,nl ! cld - do il=1,ncum ! cld - if (wa(il,i).gt.0.0) & - siga(il,i)=mac(il,i)/wa(il,i) & - *rrd*tvp(il,i)/p(il,i)/100./delta ! cld - siga(il,i) = min(siga(il,i),1.0) ! cld - !IM cf. FH - if (iflag_clw.eq.0) then - qcondc(il,i)=siga(il,i)*clw(il,i)*(1.-ep(il,i)) & - + (1.-siga(il,i))*qcond(il,i) ! cld - else if (iflag_clw.eq.1) then - qcondc(il,i)=qcond(il,i) ! cld + do i = 1, nl + do il = 1, ncum + if (i >= icb(il) .and. i <= (inb(il) - 1) & + .and. sax(il, i) > 0.0) then + wa(il, i) = sqrt(2. * sax(il, i)) endif + enddo + enddo - enddo ! cld - enddo ! cld + do i = 1, nl + do il = 1, ncum + 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) + + if (iflag_clw == 0) then + qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) & + + (1. - siga(il, i)) * qcond(il, i) + else if (iflag_clw == 1) then + qcondc(il, i) = qcond(il, i) + endif + enddo + enddo - end SUBROUTINE cv3_yield + end SUBROUTINE cv30_yield -end module cv3_yield_m +end module cv30_yield_m