Changeset 5895 for LMDZ6/trunk


Ignore:
Timestamp:
Dec 2, 2025, 4:06:07 PM (13 hours ago)
Author:
fhourdin
Message:

Nouvelles version des panaches thermiques.
Exploratoire.

Location:
LMDZ6/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_thermcell_main.F90

    r5853 r5895  
    429429         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,&
    430430     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    431      &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     431     &    lalim,f0,zmax0,detr_star,entr_star,f_star,csc,ztva,  &
    432432     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    433433     &    ,lev_out,lunout1,igout)

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

    r5512 r5895  
    11MODULE lmdz_thermcell_plume
    22!
    3 ! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $
     3! $Id: lmdz_thermcell_plume.f90 5854 2025-10-10 14:00:56Z fhourdin $
    44!
    55CONTAINS
    … …  
    77      SUBROUTINE thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    88     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    9      &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     9     &           lalim,f0,zmax0,detr_star,entr_star,f_star,csc,ztva,  &
    1010     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    1111    &           ,lev_out,lunout1,igout)
    1212!     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    1313!--------------------------------------------------------------------------
    14 ! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam
    15 !
    1614!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    17 !   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
    18 !   thermcell_plume_6A is activate for flag_thermas_ed < 10
    19 !   thermcell_plume_5B for flag_thermas_ed < 20
    20 !   thermcell_plume for flag_thermals_ed>= 20
    21 !   Various options are controled by the flag_thermals_ed parameter
    22 !   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
    23 !   = 21 : the Jam strato-cumulus modif is not activated in detrainment
    24 !   = 29 : an other way to compute the modified buoyancy (to be tested)
    2515!--------------------------------------------------------------------------
    2616
    … …  
    4636      real,intent(in),dimension(ngrid,nlay) :: pphi
    4737      real,intent(in),dimension(ngrid,nlay) :: zpspsk
    48       real,intent(in),dimension(ngrid) :: f0
     38      real,intent(in),dimension(ngrid) :: f0,zmax0
    4939
    5040      integer,intent(out) :: lalim(ngrid)
    … …  
    6454      real,intent(out),dimension(ngrid,nlay) :: ztva_est
    6555      real,intent(out),dimension(ngrid,nlay) :: zqsatth
    66       integer,intent(out),dimension(ngrid) :: lmix(ngrid)
    67       integer,intent(out),dimension(ngrid) :: lmix_bis(ngrid)
    68       real,intent(out),dimension(ngrid) :: linter(ngrid)
    69 
    70 
    71       REAL,dimension(ngrid,nlay+1) :: wa_moy
    72       REAL,dimension(ngrid,nlay) :: entr,detr
    73       REAL,dimension(ngrid,nlay) :: ztv_est
    74       REAL,dimension(ngrid,nlay) :: zqla_est
    75       REAL,dimension(ngrid,nlay) :: zta_est
    76       REAL,dimension(ngrid) :: ztemp,zqsat
     56      integer,intent(out),dimension(ngrid) :: lmix
     57      integer,intent(out),dimension(ngrid) :: lmix_bis
     58      real,intent(out),dimension(ngrid) :: linter
     59
    7760      REAL zdw2,zdw2bis
    7861      REAL zw2modif
    … …  
    8063      REAL,dimension(ngrid,nlay) :: zeps
    8164
    82       REAL,dimension(ngrid) :: wmaxa
    83 
    84       INTEGER ig,l,k,lt,it,lm,nbpb
     65      REAL, dimension(ngrid) ::    wmaxa
     66
     67      INTEGER ig,l,k,lt,it,lm
     68      integer nbpb
     69
     70      real,dimension(ngrid,nlay) :: detr
     71      real,dimension(ngrid,nlay) :: entr
     72      real,dimension(ngrid,nlay+1) :: wa_moy
     73      real,dimension(ngrid,nlay) :: ztv_est
     74      real,dimension(ngrid) :: ztemp,zqsat
     75      real,dimension(ngrid,nlay) :: zqla_est
     76      real,dimension(ngrid,nlay) :: zta_est
    8577
    8678      real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt
    8779      real zdz,zalpha,zw2m
    8880      real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam
    89       real zdz2,zdz3,lmel,entrbis,zdzbis
    90       real,dimension(ngrid) :: d_temp
     81      real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
     82      real, dimension(ngrid) :: d_temp
    9183      real ztv1,ztv2,factinv,zinv,zlmel
    9284      real zlmelup,zlmeldwn,zlt,zltdwn,zltup
    … …  
    9991      LOGICAL,dimension(ngrid) :: active,activetmp
    10092      REAL fact_gamma,fact_gamma2,fact_epsilon2
    101 
    102 
     93      REAL coefc
    10394      REAL,dimension(ngrid,nlay) :: c2
    10495
    105       if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
    106       Zsat=.false.
    107 ! Initialisation
    108 
     96      real, dimension(ngrid,nlay) :: z_detr_q,z_cld_radius
     97      real, dimension(ngrid) :: z_cld_radius_base,z_cld_base_height,z_alpha_base
     98
     99      integer choice_ed_main,choice_ed_dq
     100
     101      z_cld_base_height(:)=0.
     102
     103      if (ngrid==1) print*,'THERMCELL PLUME Modifie 2025/11/11'
     104
     105! ---------------------------------------------------------------
     106!  Controling entrainment and detrainment
     107!     iflag_thermals_ed=1XY
     108!     choice_ed_main=X
     109!     choice_ed_dq=Y
     110!     X=0 et Y=0 <=> thermcell_plume_6A
     111      choice_ed_main=(iflag_thermals_ed-100)/10
     112      choice_ed_dq=(iflag_thermals_ed-100)-10*choice_ed_main
     113! ---------------------------------------------------------------
     114 
    109115
    110116      zbetalpha=betalpha/(1.+betalpha)
    … …  
    112118
    113119! Initialisations des variables r?elles
     120Zsat=.false.
    114121if (1==1) then
    115122      ztva(:,:)=ztv(:,:)
    … …  
    154161      wmaxa(:)=0.
    155162
     163! Initialisation a 0  en cas de sortie dans replay
     164      zqsat(:)=0.
     165      zta_est(:,:)=0.
     166      zdqt(:,:)=0.
     167      zdqtjam(:,:)=0.
     168      c2(:,:)=0.
     169
    156170
    157171!-------------------------------------------------------------------------
    … …  
    221235    do ig=1,ngrid
    222236!       print*,'active',active(ig),ig,l
    223         if(active(ig)) then
     237       if(active(ig)) then
    224238        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    225239        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    226240        zta_est(ig,l)=ztva_est(ig,l)
    227241        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    228         ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    229      &      -zqla_est(ig,l))-zqla_est(ig,l))
     242        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)-zqla_est(ig,l)))
    230243 
    231244
    … …  
    250263!=========================================================================
    251264
    252 !--------------------------------------------------
    253265        lt=l+1
    254266        zlt=zlev(ig,lt)
    … …  
    256268
    257269        do while (lmel.gt.zdz2)
    258            lt=lt+1
    259            zlt=zlev(ig,lt)
    260            zdz2=zlt-zlev(ig,l)
     270          lt=lt+1
     271          zlt=zlev(ig,lt)
     272          zdz2=zlt-zlev(ig,l)
    261273        enddo
    262274        zdz3=zlev(ig,lt+1)-zlt
    … …  
    274286        zdzbis=zlev(ig,l)-zlev(ig,l-1)
    275287        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     288
    276289        zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    277290        zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    … …  
    280293        lm=Max(1,l-2)
    281294        w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
     295
     296
    282297       endif
    283298    enddo
    … …  
    285300
    286301!-------------------------------------------------
    287 !calcul des taux d'entrainement et de detrainement
     302! 4. calcul des taux d'entrainement et de detrainement
    288303!-------------------------------------------------
    289304
    … …  
    293308!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
    294309          zw2m=w_est(ig,l+1)
     310!          zw2m=zw2(ig,l)
    295311          zdz=zlev(ig,l+1)-zlev(ig,l)
    296312          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     313!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
     314          zbuoybis=zbuoy(ig,l)
    297315          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    298316          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
    299317
    300 !=========================================================================
    301 ! 4. Calcul de l'entrainement et du detrainement
    302 !=========================================================================
    303 
    304           detr_star(ig,l)=f_star(ig,l)*zdz             &
    305     &     *( mix0 * 0.1 / (zalpha+0.001)               &
    306     &     + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
    307     &     + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
    308 
    309           if ( iflag_thermals_ed == 20 ) then
    310              entr_star(ig,l)=f_star(ig,l)*zdz* (         &
    311     &          mix0 * 0.1 / (zalpha+0.001)               &
    312     &        + zbetalpha*MAX(entr_min,                   &
    313     &        afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
     318         
     319         
     320! detr_q_coef = thermals_detr_q_coef lu dans les .def
     321! vrai flux de masse : f0*fstar
     322! equation de conservation de la masse s'écrit : f*(k+1) - f*(k) = e*(k) - d*(k)
     323! e=f0 e* / dz
     324
     325          if ( choice_ed_dq == 0 )  then
     326
     327              z_detr_q(ig,l)=detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power
     328
     329          elseif ( choice_ed_dq == 1 )  then
     330
     331              ! Tant que la couche du dessous n'est pas condensée, on ne tient pas
     332              ! compte du détrainement en q
     333              ! FH V1 : if ( zqla(ig,l-1) > 0. ) then
     334              if ( zqla_est(ig,l) > 0. ) then
     335                  if ( z_cld_base_height(ig) == 0. ) then
     336                       ! Cloud base : height and fraction
     337                       z_cld_base_height(ig)=zlev(ig,l) ! z at cloud base
     338                       z_alpha_base(ig)=zalpha
     339                  endif
     340                  ! Cloud radius. At cloud base = cloud_height / 2. Cloud_height=zmax0-z_cld_base_height
     341                  !    With min value dz of the layer
     342                  ! FH V1 : z_cld_radius(ig,l)=sqrt(zalpha/z_alpha_base(ig))*0.5*(max(zmax0(ig)-z_cld_base_height(ig),zlev(ig,l+1)-zlev(ig,l)))
     343                  z_cld_radius(ig,l)=0.5*(max(zmax0(ig)-z_cld_base_height(ig),zlev(ig,l+1)-zlev(ig,l)))
     344                  z_detr_q(ig,l)=detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power*2./z_cld_radius(ig,l)
     345              else
     346                  z_detr_q(ig,l)=0.
     347              endif
     348
    314349          else
    315              entr_star(ig,l)=f_star(ig,l)*zdz* (         &
    316     &          mix0 * 0.1 / (zalpha+0.001)               &
    317     &        + zbetalpha*MAX(entr_min,                   &
    318     &        afact*zbuoy(ig,l)/zw2m - fact_epsilon))
     350
     351              print*,'choice_ed_dq=',choice_ed_dq,' non prevu'
     352              stop
     353
    319354          endif
     355
     356! zbetalpha = B1 / ( 1 + B1 )
     357! afact     = A1
     358! zbuoyjam  = B
     359! zw2m      = w2
     360! z_detr_q  = terme de détrainement de Delta q
     361
     362          if ( choice_ed_main == 0 ) then
     363
     364              detr_star(ig,l)=f_star(ig,l)*zdz             &
     365    &         *( mix0 * 0.1 / (zalpha+0.001)               &
     366    &         + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
     367    &         + z_detr_q(ig,l)) )
     368
     369              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     370              entr_star(ig,l)=f_star(ig,l)*zdz* (         &
     371    &           mix0 * 0.1 / (zalpha+0.001)               &
     372    &         + zbetalpha*MAX(entr_min,                   &
     373    &         afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
     374
     375          elseif ( choice_ed_main == 1 ) then
     376
     377              detr_star(ig,l)=f_star(ig,l)*zdz                 &
     378    &         *( mix0 * 0.1 / (zalpha+0.001)                   &
     379    &         + detr_min                                       &
     380    &         + MAX(-afact*zbetalpha*zbuoyjam(ig,l)/zw2m,0.)   &
     381    &         + z_detr_q(ig,l) )
     382
     383              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     384              entr_star(ig,l)=f_star(ig,l)*zdz* (              &
     385    &           mix0 * 0.1 / (zalpha+0.001)                    &
     386    &         + entr_min                                       &
     387    &         + zbetalpha*MAX(0.,                              &
     388    &         afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
     389
     390
     391          else
     392              print*,'choice_ed_dq=',choice_ed_dq,' non prevu'
     393              stop
     394          endif
     395
     396
    320397         
    321398! En dessous de lalim, on prend le max de alim_star et entr_star pour
    … …  
    325402          entr_star(ig,l)=0.
    326403        endif
    327         f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
     404!        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
     405!          alim_star(ig,l)=entrbis
     406!        endif
     407
     408!        print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
     409! Calcul du flux montant normalise
     410      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    328411     &              -detr_star(ig,l)
    329412
    … …  
    361444!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
    362445           zha(ig,l) = ztva(ig,l)
    363            ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
    364      &              -zqla(ig,l))-zqla(ig,l))
     446           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)-zqla(ig,l)))
    365447           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    366448           zdz=zlev(ig,l+1)-zlev(ig,l)
    … …  
    375457   enddo
    376458
    377    if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
     459        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
    378460!
    379461!===========================================================================
    … …  
    384466   do ig=1,ngrid
    385467            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    386 !               stop'On tombe sur le cas particulier de thermcell_dry'
    387 !               print*,'On tombe sur le cas particulier de thermcell_plume'
    388468                nbpb=nbpb+1
    389469                zw2(ig,l+1)=0.
    … …  
    438518        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    439519
    440 
    441520 RETURN
    442521     END SUBROUTINE thermcell_plume
Note: See TracChangeset for help on using the changeset viewer.