! $Id: lmdz_cv30.f90 5159 2024-08-02 19:58:25Z evignon $ MODULE lmdz_cv30 !------------------------------------------------------------ ! Parameters for convectL, iflag_con=30: ! (includes - microphysical parameters, ! - parameters that control the rate of approach ! to quasi-equilibrium) ! - noff & minorig (previously in input of convect1) !------------------------------------------------------------ IMPLICIT NONE; PRIVATE PUBLIC sigd, spfac, pbcrit, ptcrit, omtrain, dtovsh, dpbase, dttrig, dtcrit, & tau, beta, alpha, delta, betad, noff, minorig, nl, nlp, nlm, & cv30_param, cv30_prelim, cv30_feed, cv30_undilute1, cv30_trigger, & cv30_compress, cv30_undilute2, cv30_closure, cv30_mixing, cv30_unsat, & cv30_yield, cv30_tracer, cv30_uncompress, cv30_epmax_fn_cape INTEGER noff, minorig, nl, nlp, nlm REAL sigd, spfac REAL pbcrit, ptcrit REAL omtrain REAL dtovsh, dpbase, dttrig REAL dtcrit, tau, beta, alpha REAL delta REAL betad !$OMP THREADPRIVATE(sigd, spfac, pbcrit, ptcrit, omtrain, dtovsh, dpbase, dttrig, dtcrit, & !$OMP tau, beta, alpha, delta, betad, noff, minorig, nl, nlp, nlm) CONTAINS SUBROUTINE cv30_param(nd, delt) USE lmdz_conema3 IMPLICIT NONE ! ------------------------------------------------------------ ! Set parameters for convectL for iflag_con = 3 ! ------------------------------------------------------------ ! *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *** ! *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *** ! *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *** ! *** EFFICIENCY IS ASSUMED TO BE UNITY *** ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** ! *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** ! *** OF CLOUD *** ! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA] ! *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *** ! *** APPROACH TO QUASI-EQUILIBRIUM *** ! *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *** ! *** (BETA MUST BE LESS THAN OR EQUAL TO 1) *** ! *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *** ! *** APPROACH TO QUASI-EQUILIBRIUM *** ! *** IT MUST BE LESS THAN 0 *** INTEGER nd REAL delt ! timestep (seconds) ! noff: integer limit for convection (nd-noff) ! minorig: First level of convection ! -- limit levels for convection: noff = 1 minorig = 1 nl = nd - noff nlp = nl + 1 nlm = nl - 1 ! -- "microphysical" parameters: sigd = 0.01 spfac = 0.15 pbcrit = 150.0 ptcrit = 500.0 ! IM cf. FH epmax = 0.993 omtrain = 45.0 ! used also for snow (no disctinction rain/snow) ! -- misc: dtovsh = -0.2 ! dT for overshoot dpbase = -40. ! definition cloud base (400m above LCL) dttrig = 5. ! (loose) condition for triggering ! -- rate of approach to quasi-equilibrium: dtcrit = -2.0 tau = 8000. beta = 1.0 - delt / tau alpha = 1.5E-3 * delt / tau ! increase alpha to compensate W decrease: alpha = alpha * 1.5 ! -- interface cloud parameterization: delta = 0.01 ! cld ! -- interface with boundary-layer (gust factor): (sb) betad = 10.0 ! original value (from convect 4.3) END SUBROUTINE cv30_param SUBROUTINE cv30_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm, & th) USE lmdz_cvthermo IMPLICIT NONE ! ===================================================================== ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY ! "ori": from convect4.3 (vectorized) ! "convect3": to be exactly consistent with convect3 ! ===================================================================== ! inputs: INTEGER len, nd, ndp1 REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1) ! outputs: REAL lv(len, nd), cpn(len, nd), tv(len, nd) REAL gz(len, nd), h(len, nd), hm(len, nd) REAL th(len, nd) ! local variables: INTEGER k, i REAL rdcp REAL tvx, tvy ! convect3 REAL cpx(len, nd) ! ori do 110 k=1,nlp DO k = 1, nl ! convect3 DO i = 1, len ! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0) lv(i, k) = lv0 - clmcpv * (t(i, k) - 273.15) cpn(i, k) = cpd * (1.0 - q(i, k)) + cpv * q(i, k) cpx(i, k) = cpd * (1.0 - q(i, k)) + cl * q(i, k) ! ori tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1) tv(i, k) = t(i, k) * (1.0 + q(i, k) / eps - q(i, k)) rdcp = (rrd * (1. - q(i, k)) + q(i, k) * rrv) / cpn(i, k) th(i, k) = t(i, k) * (1000.0 / p(i, k))**rdcp END DO END DO ! gz = phi at the full levels (same as p). DO i = 1, len gz(i, 1) = 0.0 END DO ! ori do 140 k=2,nlp DO k = 2, nl ! convect3 DO i = 1, len tvx = t(i, k) * (1. + q(i, k) / eps - q(i, k)) !convect3 tvy = t(i, k - 1) * (1. + q(i, k - 1) / eps - q(i, k - 1)) !convect3 gz(i, k) = gz(i, k - 1) + 0.5 * rrd * (tvx + tvy) & !convect3 * (p(i, k - 1) - p(i, k)) / ph(i, k) !convect3 ! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k)) ! ori & *(p(i,k-1)-p(i,k))/ph(i,k) END DO END DO ! h = phi + cpT (dry static energy). ! hm = phi + cp(T-Tbase)+Lq ! ori do 170 k=1,nlp DO k = 1, nl ! convect3 DO i = 1, len h(i, k) = gz(i, k) + cpn(i, k) * t(i, k) hm(i, k) = gz(i, k) + cpx(i, k) * (t(i, k) - t(i, 1)) + lv(i, k) * q(i, k) END DO END DO END SUBROUTINE cv30_prelim SUBROUTINE cv30_feed(len, nd, t, q, qs, p, ph, hm, gz, nk, icb, icbmax, & iflag, tnk, qnk, gznk, plcl) IMPLICIT NONE ! ================================================================ ! Purpose: CONVECTIVE FEED ! Main differences with cv_feed: ! - ph added in input ! - here, nk(i)=minorig ! - icb defined differently (plcl compared with ph instead of p) ! Main differences with convect3: ! - we do not compute dplcldt and dplcldr of CLIFT anymore ! - values iflag different (but tests identical) ! - A,B explicitely defined (!...) ! ================================================================ ! inputs: INTEGER len, nd REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd) REAL hm(len, nd), gz(len, nd) REAL ph(len, nd + 1) ! outputs: INTEGER iflag(len), nk(len), icb(len), icbmax REAL tnk(len), qnk(len), gznk(len), plcl(len) ! local variables: INTEGER i, k INTEGER ihmin(len) REAL work(len) REAL pnk(len), qsnk(len), rh(len), chi(len) REAL a, b ! convect3 ! ym plcl = 0.0 ! @ !------------------------------------------------------------------- ! @ ! --- Find level of minimum moist static energy ! @ ! --- If level of minimum moist static energy coincides with ! @ ! --- or is lower than minimum allowable parcel origin level, ! @ ! --- set iflag to 6. ! @ !------------------------------------------------------------------- ! @ ! @ do 180 i=1,len ! @ work(i)=1.0e12 ! @ ihmin(i)=nl ! @ 180 continue ! @ do 200 k=2,nlp ! @ do 190 i=1,len ! @ if((hm(i,k).lt.work(i)).AND. ! @ & (hm(i,k).lt.hm(i,k-1)))THEN ! @ work(i)=hm(i,k) ! @ ihmin(i)=k ! @ endif ! @ 190 continue ! @ 200 continue ! @ do 210 i=1,len ! @ ihmin(i)=min(ihmin(i),nlm) ! @ IF(ihmin(i).le.minorig)THEN ! @ iflag(i)=6 ! @ endif ! @ 210 continue ! @ c ! @ !------------------------------------------------------------------- ! @ ! --- Find that model level below the level of minimum moist static ! @ ! --- energy that has the maximum value of moist static energy ! @ !------------------------------------------------------------------- ! @ ! @ do 220 i=1,len ! @ work(i)=hm(i,minorig) ! @ nk(i)=minorig ! @ 220 continue ! @ do 240 k=minorig+1,nl ! @ do 230 i=1,len ! @ if((hm(i,k).gt.work(i)).AND.(k.le.ihmin(i)))THEN ! @ work(i)=hm(i,k) ! @ nk(i)=k ! @ endif ! @ 230 continue ! @ 240 continue ! ------------------------------------------------------------------- ! --- Origin level of ascending parcels for convect3: ! ------------------------------------------------------------------- DO i = 1, len nk(i) = minorig END DO ! ------------------------------------------------------------------- ! --- Check whether parcel level temperature and specific humidity ! --- are reasonable ! ------------------------------------------------------------------- DO i = 1, len IF (((t(i, nk(i))<250.0) .OR. (q(i, nk(i))<=0.0)) & ! @ & .OR.( ! p(i,ihmin(i)).lt.400.0 ! ) ) .AND. (iflag(i)==0)) iflag(i) = 7 END DO ! ------------------------------------------------------------------- ! --- Calculate lifted condensation level of air at parcel origin level ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) ! ------------------------------------------------------------------- a = 1669.0 ! convect3 b = 122.0 ! convect3 DO i = 1, len IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002 tnk(i) = t(i, nk(i)) qnk(i) = q(i, nk(i)) gznk(i) = gz(i, nk(i)) pnk(i) = p(i, nk(i)) qsnk(i) = qs(i, nk(i)) rh(i) = qnk(i) / qsnk(i) ! ori rh(i)=min(1.0,rh(i)) ! removed for convect3 ! ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) chi(i) = tnk(i) / (a - b * rh(i) - tnk(i)) ! convect3 plcl(i) = pnk(i) * (rh(i)**chi(i)) IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag & (i) = 8 END IF ! iflag=7 END DO ! ------------------------------------------------------------------- ! --- Calculate first level above lcl (=icb) ! ------------------------------------------------------------------- ! @ do 270 i=1,len ! @ icb(i)=nlm ! @ 270 continue ! @c ! @ do 290 k=minorig,nl ! @ do 280 i=1,len ! @ if((k.ge.(nk(i)+1)).AND.(p(i,k).lt.plcl(i))) ! @ & icb(i)=min(icb(i),k) ! @ 280 continue ! @ 290 continue ! @c ! @ do 300 i=1,len ! @ if((icb(i).ge.nlm).AND.(iflag(i).EQ.0))iflag(i)=9 ! @ 300 continue DO i = 1, len icb(i) = nlm END DO ! la modification consiste a comparer plcl a ph et non a p: ! icb est defini par : ph(icb)p(icb), then icbs=icb ! * the routine above computes tvp from minorig to icbs (included). ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1) ! must be known. This is the case if icbs=icb+1, but not if icbs=icb. ! * therefore, in the case icbs=icb, we compute tvp at level icb+1 ! (tvp at other levels will be computed in cv3_undilute2.F) DO i = 1, len ticb(i) = t(i, icb(i) + 1) gzicb(i) = gz(i, icb(i) + 1) qsicb(i) = qs(i, icb(i) + 1) END DO DO i = 1, len tg = ticb(i) qg = qsicb(i) ! convect3 ! debug alv=lv0-clmcpv*(ticb(i)-t0) alv = lv0 - clmcpv * (ticb(i) - 273.15) ! First iteration. ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) s = cpd * (1. - qnk(i)) + cl * qnk(i) & ! convect3 + alv * alv * qg / (rrv * ticb(i) * ticb(i)) ! convect3 s = 1. / s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gzicb(i) ! convect3 tg = tg + s * (ah0(i) - ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112 * exp(17.67 * tc / denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori endif ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg = eps * es / (p(i, icb(i) + 1) - es * (1. - eps)) ! Second iteration. ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! ori s=1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gzicb(i) ! convect3 tg = tg + s * (ah0(i) - ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112 * exp(17.67 * tc / denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori end if ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg = eps * es / (p(i, icb(i) + 1) - es * (1. - eps)) alv = lv0 - clmcpv * (ticb(i) - 273.15) ! ori c approximation here: ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) ! ori & -gz(i,icb(i))-alv*qg)/cpd ! convect3: no approximation: tp(i, icb(i) + 1) = (ah0(i) - gz(i, icb(i) + 1) - alv * qg) / (cpd + (cl - cpd) * qnk(i)) ! ori clw(i,icb(i))=qnk(i)-qg ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i))) clw(i, icb(i) + 1) = qnk(i) - qg clw(i, icb(i) + 1) = max(0.0, clw(i, icb(i) + 1)) rg = qg / (1. - qnk(i)) ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) ! convect3: (qg utilise au lieu du vrai mixing ratio rg) tvp(i, icb(i) + 1) = tp(i, icb(i) + 1) * (1. + qg / eps - qnk(i)) !whole thing END DO END SUBROUTINE cv30_undilute1 SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, & iflag, sig, w0) IMPLICIT NONE ! ------------------------------------------------------------------- ! --- TRIGGERING ! - computes the cloud base ! - triggering (crude in this version) ! - relaxation of sig and w0 when no convection ! Caution1: if no convection, we set iflag=4 ! (it used to be 0 in convect3) ! Caution2: at this stage, tvp (and thus buoy) are know up ! through icb only! ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy ! ------------------------------------------------------------------- ! input: INTEGER len, nd INTEGER icb(len) REAL plcl(len), p(len, nd) REAL th(len, nd), tv(len, nd), tvp(len, nd) ! output: REAL pbase(len), buoybase(len) ! input AND output: INTEGER iflag(len) REAL sig(len, nd), w0(len, nd) ! local variables: INTEGER i, k REAL tvpbase, tvbase, tdif, ath, ath1 ! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy DO i = 1, len pbase(i) = plcl(i) + dpbase tvpbase = tvp(i, icb(i)) * (pbase(i) - p(i, icb(i) + 1)) / & (p(i, icb(i)) - p(i, icb(i) + 1)) + tvp(i, icb(i) + 1) * (p(i, icb(i)) - pbase(i)) / (& p(i, icb(i)) - p(i, icb(i) + 1)) tvbase = tv(i, icb(i)) * (pbase(i) - p(i, icb(i) + 1)) / & (p(i, icb(i)) - p(i, icb(i) + 1)) + tv(i, icb(i) + 1) * (p(i, icb(i)) - pbase(i)) / (p & (i, icb(i)) - p(i, icb(i) + 1)) buoybase(i) = tvpbase - tvbase END DO ! *** make sure that column is dry adiabatic between the surface *** ! *** and cloud base, and that lifted air is positively buoyant *** ! *** at cloud base *** ! *** if not, return to calling program after resetting *** ! *** sig(i) and w0(i) *** ! oct3 do 200 i=1,len ! oct3 ! oct3 tdif = buoybase(i) ! oct3 ath1 = th(i,1) ! oct3 ath = th(i,icb(i)-1) - dttrig ! oct3 ! oct3 if (tdif.lt.dtcrit .OR. ath.gt.ath1) THEN ! oct3 do 60 k=1,nl ! oct3 sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif ! oct3 sig(i,k) = AMAX1(sig(i,k),0.0) ! oct3 w0(i,k) = beta*w0(i,k) ! oct3 60 continue ! oct3 iflag(i)=4 ! pour version vectorisee ! oct3c convect3 iflag(i)=0 ! oct3cccc RETURN ! oct3 endif ! oct3 ! oct3200 continue ! -- oct3: on reecrit la boucle 200 (pour la vectorisation) DO k = 1, nl DO i = 1, len tdif = buoybase(i) ath1 = th(i, 1) ath = th(i, icb(i) - 1) - dttrig IF (tdifath1) THEN sig(i, k) = beta * sig(i, k) - 2. * alpha * tdif * tdif sig(i, k) = amax1(sig(i, k), 0.0) w0(i, k) = beta * w0(i, k) iflag(i) = 4 ! pour version vectorisee ! convect3 iflag(i)=0 END IF END DO END DO ! fin oct3 -- END SUBROUTINE cv30_trigger SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, & plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, & th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, & iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, & v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0) USE lmdz_print_control, ONLY: lunout USE lmdz_abort_physic, ONLY: abort_physic IMPLICIT NONE ! inputs: INTEGER len, ncum, nd, ntra, nloc INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len) REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len) REAL pbase1(len), buoybase1(len) REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd) REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd) REAL p1(len, nd), ph1(len, nd + 1), tv1(len, nd), tp1(len, nd) REAL tvp1(len, nd), clw1(len, nd) REAL th1(len, nd) REAL sig1(len, nd), w01(len, nd) REAL tra1(len, nd, ntra) ! outputs: ! en fait, on a nloc=len pour l'instant (cf cv_driver) INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc) REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc) REAL pbase(nloc), buoybase(nloc) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd) REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd) REAL p(nloc, nd), ph(nloc, nd + 1), tv(nloc, nd), tp(nloc, nd) REAL tvp(nloc, nd), clw(nloc, nd) REAL th(nloc, nd) REAL sig(nloc, nd), w0(nloc, nd) REAL tra(nloc, nd, ntra) ! local variables: INTEGER i, k, nn, j CHARACTER (LEN = 20) :: modname = 'cv30_compress' CHARACTER (LEN = 80) :: abort_message DO k = 1, nl + 1 nn = 0 DO i = 1, len IF (iflag1(i)==0) THEN nn = nn + 1 sig(nn, k) = sig1(i, k) w0(nn, k) = w01(i, k) t(nn, k) = t1(i, k) q(nn, k) = q1(i, k) qs(nn, k) = qs1(i, k) u(nn, k) = u1(i, k) v(nn, k) = v1(i, k) gz(nn, k) = gz1(i, k) h(nn, k) = h1(i, k) lv(nn, k) = lv1(i, k) cpn(nn, k) = cpn1(i, k) p(nn, k) = p1(i, k) ph(nn, k) = ph1(i, k) tv(nn, k) = tv1(i, k) tp(nn, k) = tp1(i, k) tvp(nn, k) = tvp1(i, k) clw(nn, k) = clw1(i, k) th(nn, k) = th1(i, k) END IF END DO END DO ! do 121 j=1,ntra ! do 111 k=1,nd ! nn=0 ! do 101 i=1,len ! IF(iflag1(i).EQ.0)THEN ! nn=nn+1 ! tra(nn,k,j)=tra1(i,k,j) ! END IF ! 101 continue ! 111 continue ! 121 continue IF (nn/=ncum) THEN WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum abort_message = '' CALL abort_physic(modname, abort_message, 1) END IF nn = 0 DO i = 1, len IF (iflag1(i)==0) THEN nn = nn + 1 pbase(nn) = pbase1(i) buoybase(nn) = buoybase1(i) plcl(nn) = plcl1(i) tnk(nn) = tnk1(i) qnk(nn) = qnk1(i) gznk(nn) = gznk1(i) nk(nn) = nk1(i) icb(nn) = icb1(i) icbs(nn) = icbs1(i) iflag(nn) = iflag1(i) END IF END DO END SUBROUTINE cv30_compress SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, & q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, & ep, sigp, buoy) ! epmax_cape: ajout arguments USE lmdz_conema3 USE lmdz_cvthermo IMPLICIT NONE ! --------------------------------------------------------------------- ! Purpose: ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES ! & ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD ! & ! FIND THE LEVEL OF NEUTRAL BUOYANCY ! Main differences convect3/convect4: ! - icbs (input) is the first level above LCL (may differ from icb) ! - many minor differences in the iterations ! - condensed water not removed from tvp in convect3 ! - vertical profile of buoyancy computed here (use of buoybase) ! - the determination of inb is different ! - no inb1, ONLY inb in output ! --------------------------------------------------------------------- ! inputs: INTEGER ncum, nd, nloc INTEGER icb(nloc), icbs(nloc), nk(nloc) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd) REAL p(nloc, nd) REAL tnk(nloc), qnk(nloc), gznk(nloc) REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd) REAL pbase(nloc), buoybase(nloc), plcl(nloc) ! outputs: INTEGER inb(nloc) REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd) REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd) REAL buoy(nloc, nd) ! local variables: INTEGER i, k REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit REAL by, defrac, pden REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc) LOGICAL lcape(nloc) ! ===================================================================== ! --- SOME INITIALIZATIONS ! ===================================================================== DO k = 1, nl DO i = 1, ncum ep(i, k) = 0.0 sigp(i, k) = spfac END DO END DO ! ===================================================================== ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES ! ===================================================================== ! --- The procedure is to solve the equation. ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. ! *** Calculate certain parcel quantities, including static energy *** DO i = 1, ncum ah0(i) = (cpd * (1. - qnk(i)) + cl * qnk(i)) * tnk(i) & ! debug & ! +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i) + qnk(i) * (lv0 - clmcpv * (tnk(i) - 273.15)) + gznk(i) END DO ! *** Find lifted parcel quantities above cloud base *** DO k = minorig + 1, nl DO i = 1, ncum ! ori IF(k.ge.(icb(i)+1))THEN IF (k>=(icbs(i) + 1)) THEN ! convect3 tg = t(i, k) qg = qs(i, k) ! debug alv=lv0-clmcpv*(t(i,k)-t0) alv = lv0 - clmcpv * (t(i, k) - 273.15) ! First iteration. ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) s = cpd * (1. - qnk(i)) + cl * qnk(i) & ! convect3 + alv * alv * qg / (rrv * t(i, k) * t(i, k)) ! convect3 s = 1. / s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) ! convect3 tg = tg + s * (ah0(i) - ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112 * exp(17.67 * tc / denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori endif qg = eps * es / (p(i, k) - es * (1. - eps)) ! Second iteration. ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! ori s=1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) ahg = cpd * tg + (cl - cpd) * qnk(i) * tg + alv * qg + gz(i, k) ! convect3 tg = tg + s * (ah0(i) - ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112 * exp(17.67 * tc / denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori endif qg = eps * es / (p(i, k) - es * (1. - eps)) ! debug alv=lv0-clmcpv*(t(i,k)-t0) alv = lv0 - clmcpv * (t(i, k) - 273.15) ! PRINT*,'cpd dans convect2 ',cpd ! PRINT*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd' ! PRINT*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd ! ori c approximation here: ! ori ! tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd ! convect3: no approximation: tp(i, k) = (ah0(i) - gz(i, k) - alv * qg) / (cpd + (cl - cpd) * qnk(i)) clw(i, k) = qnk(i) - qg clw(i, k) = max(0.0, clw(i, k)) rg = qg / (1. - qnk(i)) ! ori tvp(i,k)=tp(i,k)*(1.+rg*epsi) ! convect3: (qg utilise au lieu du vrai mixing ratio rg): tvp(i, k) = tp(i, k) * (1. + qg / eps - qnk(i)) ! whole thing END IF 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) ! ===================================================================== ! ori do 320 k=minorig+1,nl DO k = 1, nl ! convect3 DO i = 1, ncum pden = ptcrit - pbcrit ep(i, k) = (plcl(i) - p(i, k) - pbcrit) / pden * epmax ep(i, k) = amax1(ep(i, k), 0.0) ep(i, k) = amin1(ep(i, k), epmax) sigp(i, k) = spfac ! ori IF(k.ge.(nk(i)+1))THEN ! ori tca=tp(i,k)-t0 ! ori IF(tca.ge.0.0)THEN ! ori elacrit=elcrit ! ori else ! ori elacrit=elcrit*(1.0-tca/tlcrit) ! ori endif ! ori elacrit=max(elacrit,0.0) ! ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8) ! ori ep(i,k)=max(ep(i,k),0.0 ) ! ori ep(i,k)=min(ep(i,k),1.0 ) ! ori sigp(i,k)=sigs ! ori endif END DO END DO ! ===================================================================== ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL ! --- VIRTUAL TEMPERATURE ! ===================================================================== ! dans convect3, tvp est calcule en une seule fois, et sans retirer ! l'eau condensee (~> reversible CAPE) ! ori do 340 k=minorig+1,nl ! ori do 330 i=1,ncum ! ori IF(k.ge.(icb(i)+1))THEN ! ori tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k)) ! oric PRINT*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)' ! oric PRINT*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k) ! ori endif ! ori 330 continue ! ori 340 continue ! ori do 350 i=1,ncum ! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd ! ori 350 continue DO i = 1, ncum ! convect3 tp(i, nlp) = tp(i, nl) ! convect3 END DO ! convect3 ! ===================================================================== ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only): ! ===================================================================== ! -- this is for convect3 only: ! first estimate of buoyancy: DO i = 1, ncum DO k = 1, nl buoy(i, k) = tvp(i, k) - tv(i, k) END DO END DO ! set buoyancy=buoybase for all levels below base ! for safety, set buoy(icb)=buoybase DO i = 1, ncum DO k = 1, nl IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i, k)>=pbase(i))) THEN buoy(i, k) = buoybase(i) END IF END DO ! IM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i) buoy(i, icb(i)) = buoybase(i) END DO ! -- end convect3 ! ===================================================================== ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S ! --- LEVEL OF NEUTRAL BUOYANCY ! ===================================================================== ! -- this is for convect3 only: DO i = 1, ncum inb(i) = nl - 1 END DO DO i = 1, ncum DO k = 1, nl - 1 IF ((k>=icb(i)) .AND. (buoy(i, k)=icb(i)) .AND. (k<=inb(i))) THEN hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k & ) END IF END DO END DO END SUBROUTINE cv30_undilute2 SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, & sig, w0, cape, m) USE lmdz_cvthermo IMPLICIT NONE ! =================================================================== ! --- CLOSURE OF CONVECT3 ! vectorization: S. Bony ! =================================================================== ! input: INTEGER ncum, nd, nloc INTEGER icb(nloc), inb(nloc) REAL pbase(nloc) REAL p(nloc, nd), ph(nloc, nd + 1) REAL tv(nloc, nd), buoy(nloc, nd) ! input/output: REAL sig(nloc, nd), w0(nloc, nd) ! output: REAL cape(nloc) REAL m(nloc, nd) ! local variables: INTEGER i, j, k, icbmax REAL deltap, fac, w, amu REAL dtmin(nloc, nd), sigold(nloc, nd) ! ------------------------------------------------------- ! -- Initialization ! ------------------------------------------------------- DO k = 1, nl DO i = 1, ncum m(i, k) = 0.0 END DO END DO ! ------------------------------------------------------- ! -- Reset sig(i) and w0(i) for i>inb and i=(inb(i) + 1))) THEN sig(i, k) = beta * sig(i, k) + 2. * alpha * buoy(i, inb(i)) * abs(buoy(i, inb(& i))) sig(i, k) = amax1(sig(i, k), 0.0) w0(i, k) = beta * w0(i, k) END IF END DO END DO ! compute icbmax: icbmax = 2 DO i = 1, ncum icbmax = max(icbmax, icb(i)) END DO ! update sig and w0 below cloud base: DO k = 1, icbmax DO i = 1, ncum IF (k<=icb(i)) THEN sig(i, k) = beta * sig(i, k) - 2. * alpha * buoy(i, icb(i)) * buoy(i, icb(i)) sig(i, k) = amax1(sig(i, k), 0.0) w0(i, k) = beta * w0(i, k) END IF END DO END DO ! IF(inb.lt.(nl-1))THEN ! do 85 i=inb+1,nl-1 ! sig(i)=beta*sig(i)+2.*alpha*buoy(inb)* ! 1 abs(buoy(inb)) ! sig(i)=amax1(sig(i),0.0) ! w0(i)=beta*w0(i) ! 85 continue ! end if ! do 87 i=1,icb ! sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb) ! sig(i)=amax1(sig(i),0.0) ! w0(i)=beta*w0(i) ! 87 continue ! ------------------------------------------------------------- ! -- Reset fractional areas of updrafts and w0 at initial time ! -- and after 10 time steps of no convection ! ------------------------------------------------------------- DO k = 1, nl - 1 DO i = 1, ncum IF (sig(i, nd)<1.5 .OR. sig(i, nd)>12.0) THEN sig(i, k) = 0.0 w0(i, k) = 0.0 END IF END DO END DO ! ------------------------------------------------------------- ! -- Calculate convective available potential energy (cape), ! -- vertical velocity (w), fractional area covered by ! -- undilute updraft (sig), and updraft mass flux (m) ! ------------------------------------------------------------- DO i = 1, ncum cape(i) = 0.0 END DO ! compute dtmin (minimum buoyancy between ICB and given level k): DO i = 1, ncum DO k = 1, nl dtmin(i, k) = 100.0 END DO END DO DO i = 1, ncum DO k = 1, nl DO j = minorig, nl IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k - & 1))) THEN dtmin(i, k) = amin1(dtmin(i, k), buoy(i, j)) END IF END DO END DO END DO ! the interval on which cape is computed starts at pbase : DO k = 1, nl DO i = 1, ncum IF ((k>=(icb(i) + 1)) .AND. (k<=inb(i))) THEN deltap = min(pbase(i), ph(i, k - 1)) - min(pbase(i), ph(i, k)) cape(i) = cape(i) + rrd * buoy(i, k - 1) * deltap / p(i, k - 1) cape(i) = amax1(0.0, cape(i)) sigold(i, k) = sig(i, k) ! dtmin(i,k)=100.0 ! do 97 j=icb(i),k-1 ! mauvaise vectorisation ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j)) ! 97 continue sig(i, k) = beta * sig(i, k) + alpha * dtmin(i, k) * abs(dtmin(i, k)) sig(i, k) = amax1(sig(i, k), 0.0) sig(i, k) = amin1(sig(i, k), 0.01) fac = amin1(((dtcrit - dtmin(i, k)) / dtcrit), 1.0) w = (1. - beta) * fac * sqrt(cape(i)) + beta * w0(i, k) amu = 0.5 * (sig(i, k) + sigold(i, k)) * w m(i, k) = amu * 0.007 * p(i, k) * (ph(i, k) - ph(i, k + 1)) / tv(i, k) w0(i, k) = w END IF END DO END DO DO i = 1, ncum w0(i, icb(i)) = 0.5 * w0(i, icb(i) + 1) m(i, icb(i)) = 0.5 * m(i, icb(i) + 1) * (ph(i, icb(i)) - ph(i, icb(i) + 1)) / & (ph(i, icb(i) + 1) - ph(i, icb(i) + 2)) sig(i, icb(i)) = sig(i, icb(i) + 1) sig(i, icb(i) - 1) = sig(i, icb(i)) END DO ! cape=0.0 ! do 98 i=icb+1,inb ! deltap = min(pbase,ph(i-1))-min(pbase,ph(i)) ! cape=cape+rrd*buoy(i-1)*deltap/p(i-1) ! dcape=rrd*buoy(i-1)*deltap/p(i-1) ! dlnp=deltap/p(i-1) ! cape=amax1(0.0,cape) ! sigold=sig(i) ! dtmin=100.0 ! do 97 j=icb,i-1 ! dtmin=amin1(dtmin,buoy(j)) ! 97 continue ! sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin) ! sig(i)=amax1(sig(i),0.0) ! sig(i)=amin1(sig(i),0.01) ! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0) ! w=(1.-beta)*fac*sqrt(cape)+beta*w0(i) ! amu=0.5*(sig(i)+sigold)*w ! m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i) ! w0(i)=w ! 98 continue ! w0(icb)=0.5*w0(icb+1) ! m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2)) ! sig(icb)=sig(icb+1) ! sig(icb-1)=sig(icb) END SUBROUTINE cv30_closure SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, & u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, & vent, sij, elij, ments, qents, traent) USE lmdz_cvthermo IMPLICIT NONE ! --------------------------------------------------------------------- ! a faire: ! - changer rr(il,1) -> qnk(il) ! - vectorisation de la partie normalisation des flux (do 789...) ! --------------------------------------------------------------------- ! inputs: INTEGER ncum, nd, na, ntra, nloc INTEGER icb(nloc), inb(nloc), nk(nloc) REAL sig(nloc, nd) REAL qnk(nloc) REAL ph(nloc, nd + 1) REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd) REAL u(nloc, nd), v(nloc, nd) REAL tra(nloc, nd, ntra) ! input of convect3 REAL lv(nloc, na), h(nloc, na), hp(nloc, na) REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na) REAL m(nloc, na) ! input of convect3 ! outputs: REAL ment(nloc, na, na), qent(nloc, na, na) REAL uent(nloc, na, na), vent(nloc, na, na) REAL sij(nloc, na, na), elij(nloc, na, na) REAL traent(nloc, nd, nd, ntra) REAL ments(nloc, nd, nd), qents(nloc, nd, nd) REAL sigij(nloc, nd, nd) ! local variables: INTEGER i, j, k, il, im, jm INTEGER num1, num2 INTEGER nent(nloc, na) REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp REAL alt, smid, sjmin, sjmax, delp, delm REAL asij(nloc), smax(nloc), scrit(nloc) REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd) REAL wgh REAL zm(nloc, na) LOGICAL lwork(nloc) ! ===================================================================== ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS ! ===================================================================== ! ori do 360 i=1,ncum*nlp DO j = 1, nl DO i = 1, ncum nent(i, j) = 0 ! in convect3, m is computed in cv3_closure ! ori m(i,1)=0.0 END DO END DO ! ori do 400 k=1,nlp ! ori do 390 j=1,nlp DO j = 1, nl DO k = 1, nl DO i = 1, ncum qent(i, k, j) = rr(i, j) uent(i, k, j) = u(i, j) vent(i, k, j) = v(i, j) elij(i, k, j) = 0.0 ! ym ment(i,k,j)=0.0 ! ym sij(i,k,j)=0.0 END DO END DO END DO ! ym ment(1:ncum, 1:nd, 1:nd) = 0.0 sij(1:ncum, 1:nd, 1:nd) = 0.0 ! do k=1,ntra ! do j=1,nd ! instead nlp ! do i=1,nd ! instead nlp ! do il=1,ncum ! traent(il,i,j,k)=tra(il,j,k) ! enddo ! enddo ! enddo ! enddo zm(:, :) = 0. ! ===================================================================== ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING ! --- FRACTION (sij) ! ===================================================================== DO i = minorig + 1, nl DO j = minorig, nl DO il = 1, ncum IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il) - & 1)) .AND. (j<=inb(il))) THEN rti = rr(il, 1) - ep(il, i) * clw(il, i) bf2 = 1. + lv(il, j) * lv(il, j) * rs(il, j) / (rrv * t(il, j) * t(il, j) * cpd) anum = h(il, j) - hp(il, i) + (cpv - cpd) * t(il, j) * (rti - rr(il, j)) denom = h(il, i) - hp(il, i) + (cpd - cpv) * (rr(il, i) - rti) * t(il, j) dei = denom IF (abs(dei)<0.01) dei = 0.01 sij(il, i, j) = anum / dei sij(il, i, i) = 1.0 altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) * rti - rs(il, j) altem = altem / bf2 cwat = clw(il, j) * (1. - ep(il, j)) stemp = sij(il, i, j) IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN anum = anum - lv(il, j) * (rti - rs(il, j) - cwat * bf2) denom = denom + lv(il, j) * (rr(il, i) - rti) IF (abs(denom)<0.01) denom = 0.01 sij(il, i, j) = anum / denom altem = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) * rti - & rs(il, j) altem = altem - (bf2 - 1.) * cwat END IF IF (sij(il, i, j)>0.0 .AND. sij(il, i, j)<0.95) THEN qent(il, i, j) = sij(il, i, j) * rr(il, i) + (1. - sij(il, i, j)) * rti uent(il, i, j) = sij(il, i, j) * u(il, i) + & (1. - sij(il, i, j)) * u(il, nk(il)) vent(il, i, j) = sij(il, i, j) * v(il, i) + & (1. - sij(il, i, j)) * v(il, nk(il)) ! !!! do k=1,ntra ! !!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) ! !!! : +(1.-sij(il,i,j))*tra(il,nk(il),k) ! !!! END DO elij(il, i, j) = altem elij(il, i, j) = amax1(0.0, elij(il, i, j)) ment(il, i, j) = m(il, i) / (1. - sij(il, i, j)) nent(il, i) = nent(il, i) + 1 END IF sij(il, i, j) = amax1(0.0, sij(il, i, j)) sij(il, i, j) = amin1(1.0, sij(il, i, j)) END IF ! new END DO END DO ! do k=1,ntra ! do j=minorig,nl ! do il=1,ncum ! IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. ! : (j.ge.(icb(il)-1)).AND.(j.le.inb(il)))THEN ! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) ! : +(1.-sij(il,i,j))*tra(il,nk(il),k) ! END IF ! enddo ! enddo ! enddo ! *** if no air can entrain at level i assume that updraft detrains ! *** ! *** at that level and calculate detrained air flux and properties ! *** ! @ do 170 i=icb(il),inb(il) DO il = 1, ncum IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il, i)==0)) THEN ! @ IF(nent(il,i).EQ.0)THEN ment(il, i, i) = m(il, i) qent(il, i, i) = rr(il, nk(il)) - ep(il, i) * clw(il, i) uent(il, i, i) = u(il, nk(il)) vent(il, i, i) = v(il, nk(il)) elij(il, i, i) = clw(il, i) ! MAF sij(il,i,i)=1.0 sij(il, i, i) = 0.0 END IF END DO END DO ! do j=1,ntra ! do i=minorig+1,nl ! do il=1,ncum ! if (i.ge.icb(il) .AND. i.le.inb(il) .AND. nent(il,i).EQ.0) THEN ! traent(il,i,i,j)=tra(il,nk(il),j) ! END IF ! enddo ! enddo ! enddo DO j = minorig, nl DO i = minorig, nl DO il = 1, ncum IF ((j>=(icb(il) - 1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= & inb(il))) THEN sigij(il, i, j) = sij(il, i, j) END IF END DO END DO END DO ! @ enddo ! @170 continue ! ===================================================================== ! --- NORMALIZE ENTRAINED AIR MASS FLUXES ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING ! ===================================================================== ! ym CALL zilch(asum,ncum*nd) ! ym CALL zilch(bsum,ncum*nd) ! ym CALL zilch(csum,ncum*nd) CALL zilch(asum, nloc * nd) CALL zilch(csum, nloc * nd) CALL zilch(csum, nloc * nd) DO il = 1, ncum lwork(il) = .FALSE. END DO DO i = minorig + 1, nl num1 = 0 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1 END DO IF (num1<=0) GO TO 789 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il)) THEN lwork(il) = (nent(il, i)/=0) qp = rr(il, 1) - ep(il, i) * clw(il, i) anum = h(il, i) - hp(il, i) - lv(il, i) * (qp - rs(il, i)) + & (cpv - cpd) * t(il, i) * (qp - rr(il, i)) denom = h(il, i) - hp(il, i) + lv(il, i) * (rr(il, i) - qp) + & (cpd - cpv) * t(il, i) * (rr(il, i) - qp) IF (abs(denom)<0.01) denom = 0.01 scrit(il) = anum / denom alt = qp - rs(il, i) + scrit(il) * (rr(il, i) - qp) IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0 smax(il) = 0.0 asij(il) = 0.0 END IF END DO DO j = nl, minorig, -1 num2 = 0 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb(& il) - 1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1 END DO IF (num2<=0) GO TO 175 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb(& il) - 1) .AND. j<=inb(il) .AND. lwork(il)) THEN IF (sij(il, i, j)>1.0E-16 .AND. sij(il, i, j)<0.95) THEN wgh = 1.0 IF (j>i) THEN sjmax = amax1(sij(il, i, j + 1), smax(il)) sjmax = amin1(sjmax, scrit(il)) smax(il) = amax1(sij(il, i, j), smax(il)) sjmin = amax1(sij(il, i, j - 1), smax(il)) sjmin = amin1(sjmin, scrit(il)) IF (sij(il, i, j)<(smax(il) - 1.0E-16)) wgh = 0.0 smid = amin1(sij(il, i, j), scrit(il)) ELSE sjmax = amax1(sij(il, i, j + 1), scrit(il)) smid = amax1(sij(il, i, j), scrit(il)) sjmin = 0.0 IF (j>1) sjmin = sij(il, i, j - 1) sjmin = amax1(sjmin, scrit(il)) END IF delp = abs(sjmax - smid) delm = abs(sjmin - smid) asij(il) = asij(il) + wgh * (delp + delm) ment(il, i, j) = ment(il, i, j) * (delp + delm) * wgh END IF END IF END DO 175 END DO DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN asij(il) = amax1(1.0E-16, asij(il)) asij(il) = 1.0 / asij(il) asum(il, i) = 0.0 bsum(il, i) = 0.0 csum(il, i) = 0.0 END IF END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(& il) - 1) .AND. j<=inb(il)) THEN ment(il, i, j) = ment(il, i, j) * asij(il) END IF END DO END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(& il) - 1) .AND. j<=inb(il)) THEN asum(il, i) = asum(il, i) + ment(il, i, j) ment(il, i, j) = ment(il, i, j) * sig(il, j) bsum(il, i) = bsum(il, i) + ment(il, i, j) END IF END DO END DO DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN bsum(il, i) = amax1(bsum(il, i), 1.0E-16) bsum(il, i) = 1.0 / bsum(il, i) END IF END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(& il) - 1) .AND. j<=inb(il)) THEN ment(il, i, j) = ment(il, i, j) * asum(il, i) * bsum(il, i) END IF END DO END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb(& il) - 1) .AND. j<=inb(il)) THEN csum(il, i) = csum(il, i) + ment(il, i, j) END IF END DO END DO DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. & csum(il, i)>> DO i = 1, nd DO il = 1, ncum wdtraina(il, i) = 0.0 wdtrainm(il, i) = 0.0 END DO END DO ! RomP <<< ! *** check whether ep(inb)=0, if so, skip precipitating *** ! *** downdraft calculation *** DO il = 1, ncum lwork(il) = .TRUE. IF (ep(il, inb(il))<0.0001) lwork(il) = .FALSE. END DO CALL zilch(wdtrain, ncum) DO i = nl + 1, 1, -1 num1 = 0 DO il = 1, ncum IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1 END DO IF (num1<=0) GO TO 400 ! *** integrate liquid water equation to find condensed water *** ! *** and condensed water flux *** ! *** begin downdraft loop *** ! *** calculate detrained precipitation *** 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) wdtraina(il, i) = wdtrain(il) / grav ! Pa 26/08/10 RomP ELSE wdtrain(il) = 10.0 * ep(il, i) * m(il, i) * clw(il, i) wdtraina(il, i) = wdtrain(il) / 10. ! Pa 26/08/10 RomP END IF END IF END DO IF (i>1) THEN DO j = 1, i - 1 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.0) IF (cvflag_grav) THEN wdtrain(il) = wdtrain(il) + grav * awat * ment(il, j, i) ELSE wdtrain(il) = wdtrain(il) + 10.0 * awat * ment(il, j, i) END IF END IF END DO END DO DO il = 1, ncum IF (cvflag_grav) THEN wdtrainm(il, i) = wdtrain(il) / grav - wdtraina(il, i) ! Pm 26/08/10 RomP ELSE wdtrainm(il, i) = wdtrain(il) / 10. - wdtraina(il, i) ! Pm 26/08/10 RomP END IF END DO END IF ! *** find rain water and evaporation using provisional *** ! *** estimates of rp(i)and rp(i-1) *** DO il = 1, ncum IF (i<=inb(il) .AND. lwork(il)) THEN wt(il, i) = 45.0 IF (iph(i), pr1=1 & pr2=0 ! pour ph(i+1)0.0) THEN revap = 0.5 * (-b6 + sqrt(b6 * b6 + 4. * c6)) evap(il, i) = sigt * afac * revap water(il, i) = revap * revap ELSE evap(il, i) = -evap(il, i + 1) + 0.02 * (wdtrain(il) + sigd * wt(il, i) * & water(il, i + 1)) / (sigd * (ph(il, i) - ph(il, i + 1))) END IF ! *** calculate precipitating downdraft mass flux under *** ! *** hydrostatic approximation *** IF (i/=1) THEN tevap = amax1(0.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 END IF ! *** if hydrostatic assumption fails, *** ! *** solve cubic difference equation for downdraft theta *** ! *** and mass flux from two simultaneous differential eqns *** amfac = sigd * sigd * 70.0 * 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 xf = 100.0 * sigd * sigd * sigd * (ph(il, i) - ph(il, i + 1)) tf = b(il, i) - 5.0 * (th(il, i) - th(il, i - 1)) * t(il, i) / (lvcp(il, i) * & sigd * th(il, i)) af = xf * tf + mp(il, i + 1) * mp(il, i + 1) * tinv bf = 2. * (tinv * mp(il, i + 1))**3 + tinv * mp(il, i + 1) * xf * tf + & 50. * (p(il, i - 1) - p(il, i)) * xf * tevap fac2 = 1.0 IF (bf<0.0) fac2 = -1.0 bf = abs(bf) ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv IF (ur>=0.0) THEN sru = sqrt(ur) fac = 1.0 IF ((0.5 * bf - sru)<0.0) fac = -1.0 mp(il, i) = mp(il, i + 1) * tinv + (0.5 * bf + sru)**tinv + & fac * (abs(0.5 * bf - sru))**tinv ELSE d = atan(2. * sqrt(-ur) / (bf + 1.0E-28)) IF (fac2<0.0) d = 3.14159 - d mp(il, i) = mp(il, i + 1) * tinv + 2. * sqrt(af * tinv) * cos(d * tinv) END IF mp(il, i) = amax1(0.0, mp(il, i)) IF (cvflag_grav) THEN ! jyg : 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.0 * (p(il, i - 1) - p(il, i)) * tevap / (mp(il, & i) + sigd * 0.1) - 10.0 * (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.0 * (p(il, i - 1) - p(il, i)) * tevap / (mp(il, & i) + sigd * 0.1) - 10.0 * (th(il, i) - th(il, i - 1)) * t(il, i) / (lvcp(il, i & ) * sigd * th(il, i)) END IF b(il, i - 1) = amax1(b(il, i - 1), 0.0) END IF ! *** limit magnitude of mp(i) to meet cfl condition ! *** ampmax = 2.0 * (ph(il, i) - ph(il, i + 1)) * delti amp2 = 2.0 * (ph(il, i - 1) - ph(il, i)) * delti ampmax = amin1(ampmax, amp2) mp(il, i) = amin1(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))) END IF END IF ! i.EQ.1 ! *** find mixing ratio of precipitating downdraft *** IF (i/=inb(il)) THEN rp(il, i) = rr(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)) END IF rp(il, i) = rp(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) vp(il, i) = vp(il, i + 1) * mp(il, i + 1) + v(il, i) * (mp(il, i) - mp(il, i + & 1)) vp(il, i) = vp(il, i) / mp(il, i) ! do j=1,ntra ! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) ! testmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1)) ! : +tra(il,i,j)*(mp(il,i)-mp(il,i+1)) ! trap(il,i,j)=trap(il,i,j)/mp(il,i) ! END DO ELSE IF (mp(il, i + 1)>1.0E-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) END IF up(il, i) = up(il, i + 1) vp(il, i) = vp(il, i + 1) ! do j=1,ntra ! trap(il,i,j)=trap(il,i+1,j) ! END DO END IF END IF rp(il, i) = amin1(rp(il, i), rs(il, i)) rp(il, i) = amax1(rp(il, i), 0.0) END IF END IF END DO 400 END DO END SUBROUTINE cv30_unsat SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, & tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, & wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, & tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, & mike, tls, tps, qcondc, wd) USE lmdz_conema3 USE lmdz_cvflag USE lmdz_cvthermo IMPLICIT NONE ! inputs: INTEGER ncum, nd, na, ntra, nloc INTEGER icb(nloc), inb(nloc) REAL delt REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd) REAL tra(nloc, nd, ntra), 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), trap(nloc, nd, ntra) 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 traent(nloc, na, na, ntra) REAL tv(nloc, nd), tvp(nloc, nd) ! input/output: INTEGER iflag(nloc) ! outputs: REAL precip(nloc) REAL vprecip(nloc, nd + 1) REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd) REAL ftra(nloc, nd, ntra) 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 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 ! ------------------------------------------------------------- ! initialization: delti = 1.0 / delt DO il = 1, ncum precip(il) = 0.0 wd(il) = 0.0 ! gust vprecip(il, nd + 1) = 0. END DO 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 END DO END DO ! do j=1,ntra ! do i=1,nd ! do il=1,ncum ! ftra(il,i,j)=0.0 ! enddo ! enddo ! enddo DO i = 1, nl DO il = 1, ncum lvcp(il, i) = lv(il, i) / cpn(il, i) END DO END DO ! *** calculate surface precipitation in mm/day *** DO il = 1, ncum IF (ep(il, inb(il))>=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. END IF END IF END DO ! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s === ! MAF rajout pour lessivage DO k = 1, nl DO il = 1, ncum IF (k<=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. END IF END IF END DO END DO ! *** Calculate downdraft velocity scale *** ! *** NE PAS UTILISER POUR L'INSTANT *** ! do il=1,ncum ! wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il)) ! : /(sigd*p(il,icb(il))) ! enddo ! *** 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 END DO DO k = 2, nl DO il = 1, ncum IF (k<=inb(il)) THEN am(il) = am(il) + m(il, k) END IF END DO END DO DO il = 1, ncum ! convect3 if((0.1*dpinv*am).ge.delti)iflag(il)=4 IF (cvflag_grav) THEN IF ((0.01 * grav * work(il) * am(il))>=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))>=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)) END IF 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) END IF 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))) END IF ! cvflag_grav END DO ! il ! do j=1,ntra ! do il=1,ncum ! if (cvflag_grav) THEN ! ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il) ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) ! : +am(il)*(tra(il,2,j)-tra(il,1,j))) ! else ! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il) ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) ! : +am(il)*(tra(il,2,j)-tra(il,1,j))) ! END IF ! enddo ! enddo DO j = 2, nl DO il = 1, ncum IF (j<=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)) END IF ! cvflag_grav END IF ! j END DO END DO ! do k=1,ntra ! do j=2,nl ! do il=1,ncum ! if (j.le.inb(il)) THEN ! if (cvflag_grav) THEN ! ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1) ! : *(traent(il,j,1,k)-tra(il,1,k)) ! else ! ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1) ! : *(traent(il,j,1,k)-tra(il,1,k)) ! END IF ! END IF ! enddo ! 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<=inb(il)) num1 = num1 + 1 END DO IF (num1<=0) GO TO 500 CALL zilch(amp1, ncum) CALL zilch(ad, ncum) 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) END IF END DO END DO 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) END IF 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<=inb(il) .AND. j<=inb(il)) THEN ad(il) = ad(il) + ment(il, j, k) END IF END DO END DO 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) ! convect3 if((0.1*dpinv*amp1).ge.delti)iflag(il)=4 IF (cvflag_grav) THEN IF ((0.01 * grav * dpinv * amp1(il))>=delti) iflag(il) = 1 ! vecto ELSE IF ((0.1 * dpinv * amp1(il))>=delti) iflag(il) = 1 ! vecto END IF 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 END IF ! 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))) END IF ! cvflag_grav END IF ! i END DO ! do k=1,ntra ! 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) ! if (cvflag_grav) THEN ! ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) ! else ! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) ! END IF ! END IF ! enddo ! enddo 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) 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)) END IF ! 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 END IF ! i END DO END DO ! do j=1,ntra ! 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) ! if (cvflag_grav) THEN ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! else ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! END IF ! END IF ! enddo ! enddo ! enddo 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) 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)) END IF ! cvflag_grav END IF ! i and k END DO END DO ! do j=1,ntra ! 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 ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! else ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! END IF ! END IF ! i and k ! enddo ! enddo ! enddo 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 END IF ! cvflag_grav END IF ! i END DO ! sb: interface with the cloud parameterization: ! cld DO k = i + 1, nl DO il = 1, ncum IF (k<=inb(il) .AND. i<=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 END IF ! cld END DO ! cld END DO ! cld ! (particular case: no detraining level is found) ! cld DO il = 1, ncum ! cld IF (i<=inb(il) .AND. nent(il, i)==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 END IF ! cld END DO ! cld DO il = 1, ncum ! cld IF (i<=inb(il) .AND. nqcond(il, i)/=0.) THEN ! cld qcond(il, i) = qcond(il, i) / nqcond(il, i) ! cld END IF ! cld END DO ! do j=1,ntra ! 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) ! if (cvflag_grav) THEN ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) ! else ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) ! END IF ! END IF ! i ! enddo ! enddo 500 END DO ! *** 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))) END DO ! do j=1,ntra ! do il=1,ncum ! ex=0.1*ment(il,inb(il),inb(il)) ! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) ! : /(ph(il,inb(il))-ph(il,inb(il)+1)) ! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex ! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j) ! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) ! : /(ph(il,inb(il)-1)-ph(il,inb(il))) ! enddo ! enddo ! *** 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 END DO 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) + (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) END IF END DO END DO ! !!! do 700 i=1,icb(il)-1 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) END IF END DO END DO ! *** reset counter and return *** DO il = 1, ncum sig(il, nd) = 2.0 END DO DO i = 1, nd DO il = 1, ncum upwd(il, i) = 0.0 dnwd(il, i) = 0.0 END DO END DO DO i = 1, nl DO il = 1, ncum dnwd0(il, i) = -mp(il, i) END DO END DO DO i = nl + 1, nd DO il = 1, ncum dnwd0(il, i) = 0. END DO END DO 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 END IF END DO END DO DO i = 1, nl DO k = 1, nl DO il = 1, ncum up1(il, k, i) = 0.0 dn1(il, k, i) = 0.0 END DO END DO END DO 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) END IF END DO END DO END DO END DO 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<=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) END IF END DO END DO END DO ! !!! DO il=1,ncum ! !!! do i=icb(il),inb(il) ! !!! ! !!! upwd(il,i)=0.0 ! !!! dnwd(il,i)=0.0 ! !!! do k=i,inb(il) ! !!! up1=0.0 ! !!! dn1=0.0 ! !!! do n=1,i-1 ! !!! up1=up1+ment(il,n,k) ! !!! dn1=dn1-ment(il,k,n) ! !!! enddo ! !!! upwd(il,i)=upwd(il,i)+m(il,k)+up1 ! !!! dnwd(il,i)=dnwd(il,i)+dn1 ! !!! enddo ! !!! enddo ! !!! ! !!! ENDDO ! 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) END DO END DO DO i = nl + 1, nd DO il = 1, ncum mike(il, i) = 0. END DO END DO DO i = 1, nd DO il = 1, ncum ma(il, i) = 0 END DO END DO DO i = 1, nl DO j = i, nl DO il = 1, ncum ma(il, i) = ma(il, i) + m(il, j) END DO END DO END DO DO i = nl + 1, nd DO il = 1, ncum ma(il, i) = 0. END DO END DO DO i = 1, nl DO il = 1, ncum IF (i<=(icb(il) - 1)) THEN ma(il, i) = 0 END IF END DO END DO ! 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) END DO END DO 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 END DO ! *** 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 END DO ! cld END DO ! cld DO i = minorig, nl ! cld DO k = i + 1, nl + 1 ! cld DO il = 1, ncum ! cld IF (i<=inb(il) .AND. k<=(inb(il) + 1)) THEN ! cld mac(il, i) = mac(il, i) + m(il, k) ! cld END IF ! cld END DO ! cld END DO ! cld END DO ! cld DO i = 1, nl ! cld DO j = 1, i ! cld DO il = 1, ncum ! cld IF (i>=icb(il) .AND. i<=(inb(il) - 1) & ! cld .AND. j>=icb(il)) THEN ! cld sax(il, i) = sax(il, i) + rrd * (tvp(il, j) - tv(il, j)) & ! cld * (ph(il, j) - ph(il, j + 1)) / p(il, j) ! cld END IF ! cld END DO ! cld END DO ! cld END DO ! cld DO i = 1, nl ! cld DO il = 1, ncum ! cld IF (i>=icb(il) .AND. i<=(inb(il) - 1) & ! cld .AND. sax(il, i)>0.0) THEN ! cld wa(il, i) = sqrt(2. * sax(il, i)) ! cld END IF ! cld END DO ! cld END DO ! cld DO i = 1, nl ! cld DO il = 1, ncum ! cld IF (wa(il, i)>0.0) & ! cld siga(il, i) = mac(il, i) / wa(il, i) & ! cld * 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==0) THEN qcondc(il, i) = siga(il, i) * clw(il, i) * (1. - ep(il, i)) & ! cld + (1. - siga(il, i)) * qcond(il, i) ! cld ELSE IF (iflag_clw==1) THEN qcondc(il, i) = qcond(il, i) ! cld END IF END DO ! cld END DO ! cld END SUBROUTINE cv30_yield !RomP >>> SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, & d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb) IMPLICIT NONE ! inputs: INTEGER ncum, nd, na, nloc, len REAL ment(nloc, na, na), sij(nloc, na, na) REAL clw(nloc, nd), elij(nloc, na, na) REAL ep(nloc, na) INTEGER icb(nloc), inb(nloc) REAL vprecip(nloc, nd + 1) ! ouputs: REAL da(nloc, na), phi(nloc, na, na) REAL phi2(nloc, na, na) REAL d1a(nloc, na), dam(nloc, na) REAL epmlmmm(nloc, na, na), eplamm(nloc, na) ! variables pour tracer dans precip de l'AA et des mel ! local variables: INTEGER i, j, k, nam1 REAL epm(nloc, na, na) nam1 = na - 1 ! Introduced because ep is not defined for j=na ! variables d'Emanuel : du second indice au troisieme ! ---> tab(i,k,j) -> de l origine k a l arrivee j ! ment, sij, elij ! variables personnelles : du troisieme au second indice ! ---> tab(i,j,k) -> de k a j ! phi, phi2 ! initialisations DO j = 1, na DO i = 1, ncum da(i, j) = 0. d1a(i, j) = 0. dam(i, j) = 0. eplamm(i, j) = 0. END DO END DO DO k = 1, na DO j = 1, na DO i = 1, ncum epm(i, j, k) = 0. epmlmmm(i, j, k) = 0. phi(i, j, k) = 0. phi2(i, j, k) = 0. END DO END DO END DO ! fraction deau condensee dans les melanges convertie en precip : epm ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz DO j = 1, nam1 DO k = 1, j - 1 DO i = 1, ncum IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN !jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j) epm(i, j, k) = 1. - (1. - ep(i, j)) * clw(i, j) / max(elij(i, k, j), 1.E-16) epm(i, j, k) = max(epm(i, j, k), 0.0) END IF END DO END DO END DO DO j = 1, nam1 DO k = 1, nam1 DO i = 1, ncum IF (k>=icb(i) .AND. k<=inb(i)) THEN eplamm(i, j) = eplamm(i, j) + ep(i, j) * clw(i, j) * ment(i, j, k) * (1. - & sij(i, j, k)) END IF END DO END DO END DO DO j = 1, nam1 DO k = 1, j - 1 DO i = 1, ncum IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN epmlmmm(i, j, k) = epm(i, j, k) * elij(i, k, j) * ment(i, k, j) END IF END DO END DO END DO ! matrices pour calculer la tendance des concentrations dans cvltr.F90 DO j = 1, nam1 DO k = 1, nam1 DO i = 1, ncum da(i, j) = da(i, j) + (1. - sij(i, k, j)) * ment(i, k, j) phi(i, j, k) = sij(i, k, j) * ment(i, k, j) d1a(i, j) = d1a(i, j) + ment(i, k, j) * ep(i, k) * (1. - sij(i, k, j)) END DO END DO END DO DO j = 1, nam1 DO k = 1, j - 1 DO i = 1, ncum dam(i, j) = dam(i, j) + ment(i, k, j) * epm(i, j, k) * (1. - ep(i, k)) * (1. - & sij(i, k, j)) phi2(i, j, k) = phi(i, j, k) * epm(i, j, k) END DO END DO END DO END SUBROUTINE cv30_tracer ! RomP <<< SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, & vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, & dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, & epmlmmm, eplamm, wdtraina, wdtrainm, epmax_diag, iflag1, precip1, vprecip1, evap1, & ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, & dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, & elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1, epmax_diag1) ! epmax_cape IMPLICIT NONE ! inputs: INTEGER len, ncum, nd, ntra, nloc INTEGER idcum(nloc) INTEGER iflag(nloc) INTEGER inb(nloc) REAL precip(nloc) REAL vprecip(nloc, nd + 1), evap(nloc, nd) REAL ep(nloc, nd) REAL sig(nloc, nd), w0(nloc, nd) REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd) REAL ftra(nloc, nd, ntra) REAL ma(nloc, nd) REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd) REAL qcondc(nloc, nd) REAL wd(nloc), cape(nloc) REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd) REAL epmax_diag(nloc) ! epmax_cape ! RomP >>> REAL phi2(nloc, nd, nd) REAL d1a(nloc, nd), dam(nloc, nd) REAL wdtraina(nloc, nd), wdtrainm(nloc, nd) REAL sij(nloc, nd, nd) REAL elij(nloc, nd, nd), clw(nloc, nd) REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd) ! RomP <<< ! outputs: INTEGER iflag1(len) INTEGER inb1(len) REAL precip1(len) REAL vprecip1(len, nd + 1), evap1(len, nd) !<<< RomP REAL ep1(len, nd) !<<< RomP REAL sig1(len, nd), w01(len, nd) REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd) REAL ftra1(len, nd, ntra) REAL ma1(len, nd) REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd) REAL qcondc1(nloc, nd) REAL wd1(nloc), cape1(nloc) REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd) REAL epmax_diag1(len) ! epmax_cape ! RomP >>> REAL phi21(len, nd, nd) REAL d1a1(len, nd), dam1(len, nd) REAL wdtraina1(len, nd), wdtrainm1(len, nd) REAL sij1(len, nd, nd) REAL elij1(len, nd, nd), clw1(len, nd) REAL epmlmmm1(len, nd, nd), eplamm1(len, nd) ! RomP <<< ! local variables: INTEGER i, k, j DO i = 1, ncum precip1(idcum(i)) = precip(i) iflag1(idcum(i)) = iflag(i) wd1(idcum(i)) = wd(i) inb1(idcum(i)) = inb(i) cape1(idcum(i)) = cape(i) epmax_diag1(idcum(i)) = epmax_diag(i) ! epmax_cape END DO DO k = 1, nl DO i = 1, ncum vprecip1(idcum(i), k) = vprecip(i, k) evap1(idcum(i), k) = evap(i, k) !<<< RomP sig1(idcum(i), k) = sig(i, k) w01(idcum(i), k) = w0(i, k) ft1(idcum(i), k) = ft(i, k) fq1(idcum(i), k) = fq(i, k) fu1(idcum(i), k) = fu(i, k) fv1(idcum(i), k) = fv(i, k) ma1(idcum(i), k) = ma(i, k) upwd1(idcum(i), k) = upwd(i, k) dnwd1(idcum(i), k) = dnwd(i, k) dnwd01(idcum(i), k) = dnwd0(i, k) qcondc1(idcum(i), k) = qcondc(i, k) da1(idcum(i), k) = da(i, k) mp1(idcum(i), k) = mp(i, k) ! RomP >>> ep1(idcum(i), k) = ep(i, k) d1a1(idcum(i), k) = d1a(i, k) dam1(idcum(i), k) = dam(i, k) clw1(idcum(i), k) = clw(i, k) eplamm1(idcum(i), k) = eplamm(i, k) wdtraina1(idcum(i), k) = wdtraina(i, k) wdtrainm1(idcum(i), k) = wdtrainm(i, k) ! RomP <<< END DO END DO DO i = 1, ncum sig1(idcum(i), nd) = sig(i, nd) END DO ! do 2100 j=1,ntra ! do 2110 k=1,nd ! oct3 ! do 2120 i=1,ncum ! ftra1(idcum(i),k,j)=ftra(i,k,j) ! 2120 continue ! 2110 continue ! 2100 continue DO j = 1, nd DO k = 1, nd DO i = 1, ncum sij1(idcum(i), k, j) = sij(i, k, j) phi1(idcum(i), k, j) = phi(i, k, j) phi21(idcum(i), k, j) = phi2(i, k, j) elij1(idcum(i), k, j) = elij(i, k, j) epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j) END DO END DO END DO END SUBROUTINE cv30_uncompress SUBROUTINE cv30_epmax_fn_cape(nloc, ncum, nd & , cape, ep, hp, icb, inb, clw, nk, t, h, lv & , epmax_diag) USE lmdz_abort_physic, ONLY: abort_physic USE lmdz_conema3 USE lmdz_cvthermo IMPLICIT NONE ! On fait varier epmax en fn de la cape ! Il faut donc recalculer ep, et hp qui a déjà été calculé et ! qui en dépend ! Toutes les autres variables fn de ep sont calculées plus bas. ! inputs: INTEGER ncum, nd, nloc INTEGER icb(nloc), inb(nloc) REAL cape(nloc) REAL clw(nloc, nd), lv(nloc, nd), t(nloc, nd), h(nloc, nd) INTEGER nk(nloc) ! inouts: REAL ep(nloc, nd) REAL hp(nloc, nd) ! outputs ou local REAL epmax_diag(nloc) ! locals INTEGER i, k REAL hp_bak(nloc, nd) CHARACTER (LEN = 20) :: modname = 'cv30_epmax_fn_cape' CHARACTER (LEN = 80) :: abort_message ! on recalcule ep et hp IF (coef_epmax_cape>1e-12) THEN DO i = 1, ncum epmax_diag(i) = epmax - coef_epmax_cape * sqrt(cape(i)) DO k = 1, nl ep(i, k) = ep(i, k) / epmax * epmax_diag(i) ep(i, k) = amax1(ep(i, k), 0.0) ep(i, k) = amin1(ep(i, k), epmax_diag(i)) enddo enddo ! On recalcule hp: DO k = 1, nl DO i = 1, ncum hp_bak(i, k) = hp(i, k) enddo enddo DO k = 1, nlp DO i = 1, ncum hp(i, k) = h(i, k) enddo enddo DO k = minorig + 1, nl DO i = 1, ncum IF((k>=icb(i)).AND.(k<=inb(i)))THEN hp(i, k) = h(i, nk(i)) + (lv(i, k) + (cpd - cpv) * t(i, k)) * ep(i, k) * clw(i, k) endif enddo enddo !do k=minorig+1,n ! WRITE(*,*) 'cv30_routines 6218: hp(1,20)=',hp(1,20) DO i = 1, ncum DO k = 1, nl IF (abs(hp_bak(i, k) - hp(i, k))>0.01) THEN WRITE(*, *) 'i,k=', i, k WRITE(*, *) 'coef_epmax_cape=', coef_epmax_cape WRITE(*, *) 'epmax_diag(i)=', epmax_diag(i) WRITE(*, *) 'ep(i,k)=', ep(i, k) WRITE(*, *) 'hp(i,k)=', hp(i, k) WRITE(*, *) 'hp_bak(i,k)=', hp_bak(i, k) WRITE(*, *) 'h(i,k)=', h(i, k) WRITE(*, *) 'nk(i)=', nk(i) WRITE(*, *) 'h(i,nk(i))=', h(i, nk(i)) WRITE(*, *) 'lv(i,k)=', lv(i, k) WRITE(*, *) 't(i,k)=', t(i, k) WRITE(*, *) 'clw(i,k)=', clw(i, k) WRITE(*, *) 'cpd,cpv=', cpd, cpv CALL abort_physic(modname, abort_message, 1) endif enddo !do k=1,nl enddo !do i=1,ncum ENDIF !if (coef_epmax_cape.gt.1e-12) THEN END SUBROUTINE cv30_epmax_fn_cape END MODULE lmdz_cv30