Index: LMDZ6/trunk/libf/phylmd/cv3p1_closure.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/cv3p1_closure.f90 (revision 6034) +++ LMDZ6/trunk/libf/phylmd/cv3p1_closure.f90 (revision 6041) @@ -165,4 +165,8 @@ !ym column independance !ym DO k = 1, icbmax + +! JYCF2026/01/20 m(k)=rho(k)*sig(k)*w(k) +! sig(k) est la fraction surfacique de l'ascendance arrivant au niveau k (en m^2/m^2) + DO k = 1, nd DO il = 1, ncum @@ -363,4 +367,10 @@ IF (ok_inhib) THEN + +! JYCF2026/01/20 m(k)=rho(k)*sig(k)*w(k) +! Possible que toutes les lignes au dessus ne servent que pour ok_inhib=T +! Ce n'est pas le cas. +! ok_inhib est mis à iflag_mix==2 dans cva_driver. +! Il est donc faux par défaut, pour iflag_mix=1 DO i = 1, nl @@ -413,11 +423,21 @@ cinb, plfc) +! JYCF2026/01/20 somme de la partie en dessous et au dessus du LCL + DO il = 1, ncum cin(il) = cina(il) + cinb(il) END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine' + ! ------------------------------------------------------------- ! --Update buoyancies to account for Ale ! ------------------------------------------------------------- + +! JYCF2026/01/20 Ajout de Ale dans le caclul de ALE pour le déclenchement +! Est-ce que le calcul du déclenchement est fait dedans. +! A verifier et documenter. +! Notamment mettre les inten(in/out) dans cv3_buoy.f90 +! On ne sort buoy que si ALE+CIN>0 ? +! Le "buoy" qui sort est un delta de température ? CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, & @@ -436,5 +456,5 @@ ! compute dtmin (minimum buoyancy between ICB and given level k): - + ! in fact a delta of temperature DO k = 1, nl DO il = 1, ncum @@ -454,5 +474,10 @@ END DO - ! the interval on which cape is computed starts at pbase : + ! the vertical interval on which cape is computed starts at pbase : + +! JYCF2026/01/20 computing siglim(k), the target value for sig(k) +! for the time relaxation +! mlim=rho*dp*siglim*wlim avec wlim=sqrt(cape) + DO k = 1, nl @@ -498,4 +523,8 @@ END IF +! JYCF2026/01/20 computing mlim at CB +! Boucle à vérifier +! inb est le sommet + IF (icb(il)+1<=inb(il)) THEN ! IM end @@ -514,4 +543,6 @@ ! ------------------------------------------------------------------------ + +! JYCF2026/01/20 cbmflim=sum(mlim) DO il = 1, ncum cbmflim(il) = 0. @@ -541,4 +572,7 @@ wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.)) END DO + +! JYCF2026/01/20 Various options to compute the vertical velocity +! in the adiabatic ascent at cloud base DO il = 1, ncum @@ -572,4 +606,7 @@ ENDDO !RC + +! JYCF2026/01/20 computing cbmf1, mass flux at cloud base coming +! from the ALP closure DO il = 1, ncum @@ -592,4 +629,5 @@ END DO +! JYCF2026/01/20 cbmf=cbmf1, coming from ALP closure DO il = 1, ncum IF (cbmflim(il)>1.E-6) THEN @@ -607,4 +645,7 @@ ! c 2. Compute coefficient and apply correction +! JYCF2026/01/20 coef : ratio of the ALP to the CAPE closure +! keeping the original ratio in coeftrue + DO il = 1, ncum coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10) @@ -612,4 +653,21 @@ END DO IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS' + +! JYCF2026/01/20 : resumé +! computing siglim(k), the target value for sig(k) +! for the time relaxation +! mlim=rho*dp*siglim*wlim avec wlim=sqrt(cape) +! cbmflim=sum(mlim) +! computing cbmf, mass flux at cloud base from the ALP closure +! time relaxation on sig*w (in fact on "m") +! w is kept unchanged and sig is computed as sig*w/w +! +! coef=cbmf(ALP)/cbmflim(CAPE) sans relaxation +! + + +! JYCF2026/01/20 time relaxation on sig*w (in fact on "m") +! beta = 1.0 - delt / tau, in cv3_param +! w is kept unchanged and sig is computed as sig*w/w DO k = 1, nl @@ -623,4 +681,6 @@ sig(il, k) = min(sig(il,k), 1.) ! c amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k) +! JYCF2026/01/20 0.007=2/R ? +! m=sig w rho delta(p) m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k) END IF @@ -640,82 +700,82 @@ !computation of the sum of ascending fluxes IF (iflag_mix_adiab.eq.1) THEN - -!Verification sum(me)=sum(m) - DO k = 1,nd !jyg: initialization up to nd - DO il = 1, ncum - md(il,k)=0. - med(il,k)=0. - ENDDO - ENDDO - - DO k = nl,1,-1 - DO il = 1, ncum - md(il,k)=md(il,k+1)+m(il,k+1) - ENDDO - ENDDO - - DO k = nl,1,-1 - DO il = 1, ncum - IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN - mad(il,k)=mad(il,k+1)+m(il,k+1) - ENDIF -! print*,"mad",il,k,mad(il,k) - ENDDO - ENDDO - -!CR: erosion of each adiabatic ascent during its ascent - -!Computation of erosion coefficient beta_coef - DO k = 1, nl - DO il = 1, ncum - IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN -! print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1) - beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2 - ELSE - beta_coef(il,k)=0. - ENDIF - ENDDO - ENDDO - -! print*,"apres beta_coef" - - DO k = 1, nl - DO il = 1, ncum - - IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN - -! print*,"dz",il,k,tv(il, k-1) - dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG) - betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz) -! betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz) -! print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k) - dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG) -! me(il,k)=betalim(il,k)*(m(il,k)+RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz*mad(il,k)) - me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k)) -! print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz - - END IF - -!Modification of m - m(il,k)=me(il,k) - END DO - END DO - -! DO il = 1, ncum -! dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG) -! m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) & -! /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il)) -! print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il)) -! ENDDO - -!Verification sum(me)=sum(m) - DO k = nl,1,-1 - DO il = 1, ncum - med(il,k)=med(il,k+1)+m(il,k+1) -! print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k) - ENDDO - ENDDO - - + ! + !Verification sum(me)=sum(m) + DO k = 1,nd !jyg: initialization up to nd + DO il = 1, ncum + md(il,k)=0. + med(il,k)=0. + ENDDO + ENDDO + ! + DO k = nl,1,-1 + DO il = 1, ncum + md(il,k)=md(il,k+1)+m(il,k+1) + ENDDO + ENDDO + ! + DO k = nl,1,-1 + DO il = 1, ncum + IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN + mad(il,k)=mad(il,k+1)+m(il,k+1) + ENDIF + ! print*,"mad",il,k,mad(il,k) + ENDDO + ENDDO + ! + !CR: erosion of each adiabatic ascent during its ascent + ! + !Computation of erosion coefficient beta_coef + DO k = 1, nl + DO il = 1, ncum + IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN + ! print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1) + beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2 + ELSE + beta_coef(il,k)=0. + ENDIF + ENDDO + ENDDO + ! + ! print*,"apres beta_coef" + ! + DO k = 1, nl + DO il = 1, ncum + ! + IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN + ! + ! print*,"dz",il,k,tv(il, k-1) + dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG) + betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz) + ! betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz) + ! print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k) + dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG) + ! me(il,k)=betalim(il,k)*(m(il,k)+RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz*mad(il,k)) + me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k)) + ! print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz + ! + END IF + ! + !Modification of m + m(il,k)=me(il,k) + END DO + END DO + ! + ! DO il = 1, ncum + ! dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG) + ! m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) & + ! /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il)) + ! print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il)) + ! ENDDO + ! + !Verification sum(me)=sum(m) + DO k = nl,1,-1 + DO il = 1, ncum + med(il,k)=med(il,k+1)+m(il,k+1) + ! print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k) + ENDDO + ENDDO + ! + ! ENDIF !(iflag_mix_adiab) !RC @@ -759,7 +819,6 @@ ! c 4. Introduce a correcting factor for coef, in order to obtain an - ! effective - ! c sigdz larger in the present case (using cv3p1_closure) than in the - ! old + ! effective sigdz larger in the present case (using cv3p1_closure) + ! than in the old ! c closure (using cv3_closure). IF (1==0) THEN