Changeset 6041


Ignore:
Timestamp:
Jan 21, 2026, 3:48:50 AM (16 hours ago)
Author:
fhourdin
Message:

cv3p1_closure with additional comments (JY,Cat,Fredho)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/cv3p1_closure.f90

    r5708 r6041  
    165165!ym column independance
    166166!ym  DO k = 1, icbmax
     167
     168! JYCF2026/01/20 m(k)=rho(k)*sig(k)*w(k)
     169! sig(k) est la fraction surfacique de l'ascendance arrivant au niveau k (en m^2/m^2)
     170
    167171  DO k = 1, nd
    168172    DO il = 1, ncum
    … …  
    363367
    364368  IF (ok_inhib) THEN
     369
     370! JYCF2026/01/20 m(k)=rho(k)*sig(k)*w(k)
     371! Possible que toutes les lignes au dessus ne  servent que pour ok_inhib=T
     372! Ce n'est pas le cas.
     373! ok_inhib est mis à iflag_mix==2 dans cva_driver.
     374! Il est donc faux par défaut, pour iflag_mix=1
    365375
    366376    DO i = 1, nl
    … …  
    413423    cinb, plfc)
    414424
     425! JYCF2026/01/20 somme de la partie en dessous et au dessus du LCL
     426
    415427  DO il = 1, ncum
    416428    cin(il) = cina(il) + cinb(il)
    417429  END DO
    418430  IF (prt_level>=20) PRINT *, 'cv3p1_param apres cv3_cine'
     431
    419432  ! -------------------------------------------------------------
    420433  ! --Update buoyancies to account for Ale
    421434  ! -------------------------------------------------------------
     435
     436! JYCF2026/01/20 Ajout de Ale dans le caclul de ALE pour le déclenchement
     437!  Est-ce que le calcul du déclenchement est fait dedans.
     438!  A verifier et documenter.
     439!  Notamment mettre les inten(in/out) dans cv3_buoy.f90
     440!  On ne sort buoy que si ALE+CIN>0 ?
     441!  Le "buoy" qui sort est un delta de température ?
    422442
    423443  CALL cv3_buoy(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, ale, cin, tv, &
    … …  
    436456
    437457  ! compute dtmin (minimum buoyancy between ICB and given level k):
    438 
     458  ! in fact a delta of temperature
    439459  DO k = 1, nl
    440460    DO il = 1, ncum
    … …  
    454474  END DO
    455475
    456   ! the interval on which cape is computed starts at pbase :
     476  ! the vertical interval on which cape is computed starts at pbase :
     477
     478! JYCF2026/01/20 computing siglim(k), the target value for sig(k)
     479!   for the time relaxation
     480!   mlim=rho*dp*siglim*wlim avec wlim=sqrt(cape)
     481
    457482
    458483  DO k = 1, nl
    … …  
    498523    END IF
    499524
     525! JYCF2026/01/20 computing mlim at CB
     526! Boucle à vérifier
     527! inb est le sommet
     528
    500529    IF (icb(il)+1<=inb(il)) THEN
    501530      ! IM end
    … …  
    514543  ! ------------------------------------------------------------------------
    515544
     545
     546! JYCF2026/01/20 cbmflim=sum(mlim)
    516547  DO il = 1, ncum
    517548    cbmflim(il) = 0.
    … …  
    541572    wb2(il) = sqrt(2.*max(ale(il)+cin(il),0.))
    542573  END DO
     574
     575! JYCF2026/01/20 Various options to compute the vertical velocity
     576! in the adiabatic ascent at cloud base
    543577
    544578  DO il = 1, ncum
    … …  
    572606  ENDDO
    573607!RC
     608
     609! JYCF2026/01/20 computing cbmf1, mass flux at cloud base coming
     610! from the ALP closure
    574611
    575612  DO il = 1, ncum
    … …  
    592629  END DO
    593630
     631! JYCF2026/01/20 cbmf=cbmf1, coming from ALP closure
    594632  DO il = 1, ncum
    595633    IF (cbmflim(il)>1.E-6) THEN
    … …  
    607645  ! c 2. Compute coefficient and apply correction
    608646
     647! JYCF2026/01/20 coef : ratio of the ALP to the CAPE closure
     648! keeping the original ratio in coeftrue
     649
    609650  DO il = 1, ncum
    610651    coef(il) = (cbmf(il)+1.E-10)/(cbmflim(il)+1.E-10)
    … …  
    612653  END DO
    613654  IF (prt_level>=20) PRINT *, 'cv3p1_param apres coef_plantePLUS'
     655
     656! JYCF2026/01/20 : resumé
     657! computing siglim(k), the target value for sig(k)
     658!   for the time relaxation
     659!   mlim=rho*dp*siglim*wlim avec wlim=sqrt(cape)
     660! cbmflim=sum(mlim)
     661! computing cbmf, mass flux at cloud base from the ALP closure
     662! time relaxation on sig*w (in fact on "m")
     663! w is kept unchanged and sig is computed as sig*w/w
     664!
     665! coef=cbmf(ALP)/cbmflim(CAPE) sans relaxation
     666!
     667
     668
     669! JYCF2026/01/20 time relaxation on sig*w (in fact on "m")
     670! beta = 1.0 - delt / tau, in cv3_param
     671! w is kept unchanged and sig is computed as sig*w/w
    614672
    615673  DO k = 1, nl
    … …  
    623681        sig(il, k) = min(sig(il,k), 1.)
    624682        ! c         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
     683! JYCF2026/01/20 0.007=2/R ?
     684!  m=sig w rho delta(p)
    625685        m(il, k) = amu*0.007*p(il, k)*(ph(il,k)-ph(il,k+1))/tv(il, k)
    626686      END IF
    … …  
    640700!computation of the sum of ascending fluxes
    641701  IF (iflag_mix_adiab.eq.1) THEN
    642 
    643 !Verification sum(me)=sum(m)
    644   DO k = 1,nd                         !jyg: initialization up to nd
    645     DO il = 1, ncum
    646        md(il,k)=0.
    647        med(il,k)=0.
    648     ENDDO
    649   ENDDO
    650 
    651   DO k = nl,1,-1
    652     DO il = 1, ncum
    653            md(il,k)=md(il,k+1)+m(il,k+1)
    654     ENDDO
    655   ENDDO
    656 
    657   DO k = nl,1,-1
    658     DO il = 1, ncum
    659         IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
    660            mad(il,k)=mad(il,k+1)+m(il,k+1)
    661         ENDIF
    662 !        print*,"mad",il,k,mad(il,k)
    663     ENDDO
    664   ENDDO
    665 
    666 !CR: erosion of each adiabatic ascent during its ascent
    667 
    668 !Computation of erosion coefficient beta_coef
    669   DO k = 1, nl
    670     DO il = 1, ncum
    671        IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
    672 !          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
    673           beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
    674        ELSE
    675           beta_coef(il,k)=0.
    676        ENDIF
    677     ENDDO
    678   ENDDO
    679 
    680 !  print*,"apres beta_coef"
    681 
    682   DO k = 1, nl
    683     DO il = 1, ncum
    684 
    685       IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
    686 
    687 !        print*,"dz",il,k,tv(il, k-1)
    688         dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
    689         betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
    690 !        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
    691 !        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
    692         dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
    693 !        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))
    694         me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
    695 !        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz   
    696      
    697       END IF
    698        
    699 !Modification of m
    700       m(il,k)=me(il,k)
    701     END DO
    702   END DO
    703  
    704 !  DO il = 1, ncum
    705 !     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
    706 !     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
    707 !                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
    708 !     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
    709 !  ENDDO
    710 
    711 !Verification sum(me)=sum(m)
    712   DO k = nl,1,-1
    713     DO il = 1, ncum
    714            med(il,k)=med(il,k+1)+m(il,k+1)
    715 !           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
    716     ENDDO
    717   ENDDO
    718 
    719 
     702      !
     703      !Verification sum(me)=sum(m)
     704        DO k = 1,nd                         !jyg: initialization up to nd
     705          DO il = 1, ncum
     706             md(il,k)=0.
     707             med(il,k)=0.
     708          ENDDO
     709        ENDDO
     710      !
     711        DO k = nl,1,-1
     712          DO il = 1, ncum
     713                 md(il,k)=md(il,k+1)+m(il,k+1)
     714          ENDDO
     715        ENDDO
     716      !
     717        DO k = nl,1,-1
     718          DO il = 1, ncum
     719              IF ((k>=(icb(il))) .AND. (k<=inb(il))) THEN
     720                 mad(il,k)=mad(il,k+1)+m(il,k+1)
     721              ENDIF
     722      !        print*,"mad",il,k,mad(il,k)
     723          ENDDO
     724        ENDDO
     725      !
     726      !CR: erosion of each adiabatic ascent during its ascent
     727      !
     728      !Computation of erosion coefficient beta_coef
     729        DO k = 1, nl
     730          DO il = 1, ncum
     731             IF ((k>=(icb(il)+1)) .AND. (k<=inb(il)) .AND. (mlim(il,k).gt.0.)) THEN     
     732      !          print*,"beta_coef",il,k,icb(il),inb(il),buoy(il,k),tv(il,k),wlim(il,k),wlim(il,k+1)
     733                beta_coef(il,k)=RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2
     734             ELSE
     735                beta_coef(il,k)=0.
     736             ENDIF
     737          ENDDO
     738        ENDDO
     739      !
     740      !  print*,"apres beta_coef"
     741      !
     742        DO k = 1, nl
     743          DO il = 1, ncum
     744      !
     745            IF ((k>=(icb(il)+1)) .AND. (k<=inb(il))) THEN
     746      !
     747      !        print*,"dz",il,k,tv(il, k-1)
     748              dz = (ph(il,k-1)-ph(il,k))/(p(il, k-1)/(rrd*tv(il, k-1))*RG)
     749              betalim(il,k)=betalim(il,k-1)*exp(-1.*beta_coef(il,k-1)*dz)
     750      !        betalim(il,k)=betalim(il,k-1)*exp(-RG*coef_peel*buoy(il,k-1)/tv(il,k-1)/5.**2*dz)
     751      !        print*,"me",il,k,mlim(il,k),buoy(il,k),wlim(il,k),mad(il,k)
     752              dz = (ph(il,k)-ph(il,k+1))/(p(il, k)/(rrd*tv(il, k))*RG)
     753      !        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))
     754              me(il,k)=betalim(il,k)*(m(il,k)+beta_coef(il,k)*dz*mad(il,k))
     755      !        print*,"B/w2",il,k,RG*coef_peel*buoy(il,k)/tv(il,k)/((wlim(il,k)+wlim(il,k+1))/2.)**2*dz   
     756      !
     757            END IF
     758      !
     759      !Modification of m
     760            m(il,k)=me(il,k)
     761          END DO
     762        END DO
     763      !
     764      !  DO il = 1, ncum
     765      !     dz = (ph(il,icb(il))-ph(il,icb(il)+1))/(p(il, icb(il))/(rrd*tv(il, icb(il)))*RG)
     766      !     m(il,icb(il))=m(il,icb(il))+RG*coef_peel*buoy(il,icb(il))/tv(il,icb(il)) &
     767      !                  /((wlim(il,icb(il))+wlim(il,icb(il)+1))/2.)**2*dz*mad(il,icb(il))
     768      !     print*,"wlim(icb)",icb(il),wlim(il,icb(il)),m(il,icb(il))
     769      !  ENDDO
     770      !
     771      !Verification sum(me)=sum(m)
     772        DO k = nl,1,-1
     773          DO il = 1, ncum
     774                 med(il,k)=med(il,k+1)+m(il,k+1)
     775      !           print*,"somme(me),somme(m)",il,k,icb(il),med(il,k),md(il,k),me(il,k),m(il,k),wlim(il,k)
     776          ENDDO
     777        ENDDO
     778      !
     779      !
    720780  ENDIF !(iflag_mix_adiab)
    721781!RC
    … …  
    759819
    760820  ! c 4. Introduce a correcting factor for coef, in order to obtain an
    761   ! effective
    762   ! c    sigdz larger in the present case (using cv3p1_closure) than in the
    763   ! old
     821  ! effective sigdz larger in the present case (using cv3p1_closure)
     822  !  than in the old
    764823  ! c    closure (using cv3_closure).
    765824  IF (1==0) THEN
Note: See TracChangeset for help on using the changeset viewer.