Ignore:
Timestamp:
Jan 14, 2010, 3:35:30 PM (14 years ago)
Author:
Laurent Fairhead
Message:

Modifications pour la nouvelle version des thermiques (2009/2010) CR et FH

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90

    r1057 r1294  
     1!
     2! $Header$
     3!
    14      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    2      &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
    3      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    4      &           ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
     5     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
     6     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     7     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    58     &           ,lev_out,lunout1,igout)
    69
    … …  
    3134      REAL zpspsk(ngrid,klev)
    3235      REAL alim_star(ngrid,klev)
    33       REAL zmax_sec(ngrid)
    3436      REAL f0(ngrid)
    35       REAL l_mix
    36       REAL r_aspect
    3737      INTEGER lalim(ngrid)
    3838      integer lev_out                           ! niveau pour les print
    … …  
    4444      REAL ztla(ngrid,klev)
    4545      REAL zqla(ngrid,klev)
    46       REAL zqla0(ngrid,klev)
    4746      REAL zqta(ngrid,klev)
    4847      REAL zha(ngrid,klev)
    … …  
    5049      REAL detr_star(ngrid,klev)
    5150      REAL coefc
    52       REAL detr_stara(ngrid,klev)
    53       REAL detr_starb(ngrid,klev)
    54       REAL detr_starc(ngrid,klev)
    55       REAL detr_star0(ngrid,klev)
    56       REAL detr_star1(ngrid,klev)
    57       REAL detr_star2(ngrid,klev)
    58 
    5951      REAL entr_star(ngrid,klev)
    60       REAL entr_star1(ngrid,klev)
    61       REAL entr_star2(ngrid,klev)
    6252      REAL detr(ngrid,klev)
    6353      REAL entr(ngrid,klev)
     54
     55      REAL csc(ngrid,klev)
    6456
    6557      REAL zw2(ngrid,klev+1)
    … …  
    7264      REAL zqsatth(ngrid,klev)
    7365      REAL zta_est(ngrid,klev)
     66      REAL zdw2
     67      REAL zw2modif
     68      REAL zeps
    7469
    7570      REAL linter(ngrid)
    … …  
    8075      INTEGER ig,l,k
    8176
     77      real zdz,zfact,zbuoy,zalpha
    8278      real zcor,zdelta,zcvm5,qlbef
    8379      real Tbef,qsatbef
    … …  
    8682      PARAMETER (DDT0=.01)
    8783      logical Zsat
    88       REAL fact_gamma,fact_epsilon
     84      LOGICAL active(ngrid),activetmp(ngrid)
     85      REAL fact_gamma,fact_epsilon,fact_gamma2
    8986      REAL c2(ngrid,klev)
    9087
     88      REAL zw2fact,expa
    9189      Zsat=.false.
    9290! Initialisation
    … …  
    9795         fact_epsilon=1.
    9896      else if (iflag_thermals_ed==1)  then
     97! Valeurs au moment de la premiere soumission des papiers
    9998         fact_gamma=1.
    100          fact_epsilon=1.
     99         fact_epsilon=0.002
     100         fact_gamma2=0.6
     101! Suggestions des Fleurs, Septembre 2009
     102         fact_epsilon=0.015
     103!test cr
     104!         fact_epsilon=0.002
     105         fact_gamma=0.9
     106         fact_gamma2=0.7
     107
    101108      else if (iflag_thermals_ed==2)  then
    102109         fact_gamma=1.
    103110         fact_epsilon=2.
     111      else if (iflag_thermals_ed==3)  then
     112         fact_gamma=3./4.
     113         fact_epsilon=3.
    104114      endif
    105115
    106       do l=1,klev
     116      write(lunout,*)'THERM 31H '
     117
     118! Initialisations des variables reeles
     119if (1==0) then
     120      ztva(:,:)=ztv(:,:)
     121      ztva_est(:,:)=ztva(:,:)
     122      ztla(:,:)=zthl(:,:)
     123      zqta(:,:)=po(:,:)
     124      zha(:,:) = ztva(:,:)
     125else
     126      ztva(:,:)=0.
     127      ztva_est(:,:)=0.
     128      ztla(:,:)=0.
     129      zqta(:,:)=0.
     130      zha(:,:) =0.
     131endif
     132
     133      zqla_est(:,:)=0.
     134      zqsatth(:,:)=0.
     135      zqla(:,:)=0.
     136      detr_star(:,:)=0.
     137      entr_star(:,:)=0.
     138      alim_star(:,:)=0.
     139      alim_star_tot(:)=0.
     140      csc(:,:)=0.
     141      detr(:,:)=0.
     142      entr(:,:)=0.
     143      zw2(:,:)=0.
     144      w_est(:,:)=0.
     145      f_star(:,:)=0.
     146      wa_moy(:,:)=0.
     147      linter(:)=1.
     148      linter(:)=1.
     149
     150! Initialisation des variables entieres
     151      lmix(:)=1
     152      lmix_bis(:)=2
     153      wmaxa(:)=0.
     154      lalim(:)=1
     155
     156!-------------------------------------------------------------------------
     157! On ne considere comme actif que les colonnes dont les deux premieres
     158! couches sont instables.
     159!-------------------------------------------------------------------------
     160      active(:)=ztv(:,1)>ztv(:,2)
     161
     162!-------------------------------------------------------------------------
     163! Definition de l'alimentation a l'origine dans thermcell_init
     164!-------------------------------------------------------------------------
     165      do l=1,klev-1
    107166         do ig=1,ngrid
    108             zqla_est(ig,l)=0.
    109             ztva_est(ig,l)=ztva(ig,l)
    110             zqsatth(ig,l)=0.
     167            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
     168               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
     169     &                       *sqrt(zlev(ig,l+1))
     170               lalim(:)=l+1
     171               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     172            endif
    111173         enddo
    112174      enddo
    113 
    114 !CR: attention test couche alim
    115 !     do l=2,klev
    116 !     do ig=1,ngrid
    117 !        alim_star(ig,l)=0.
    118 !     enddo
    119 !     enddo
    120 !AM:initialisations du thermique
    121       do k=1,klev
    122          do ig=1,ngrid
    123             ztva(ig,k)=ztv(ig,k)
    124             ztla(ig,k)=zthl(ig,k)
    125             zqla(ig,k)=0.
    126             zqta(ig,k)=po(ig,k)
    127 !
    128             ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
    129             ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
    130             zha(ig,k) = ztva(ig,k)
    131 !
    132          enddo
    133       enddo
    134       do k=1,klev
    135         do ig=1,ngrid
    136            detr_star(ig,k)=0.
    137            entr_star(ig,k)=0.
    138 
    139            detr_stara(ig,k)=0.
    140            detr_starb(ig,k)=0.
    141            detr_starc(ig,k)=0.
    142            detr_star0(ig,k)=0.
    143            zqla0(ig,k)=0.
    144            detr_star1(ig,k)=0.
    145            detr_star2(ig,k)=0.
    146            entr_star1(ig,k)=0.
    147            entr_star2(ig,k)=0.
    148 
    149            detr(ig,k)=0.
    150            entr(ig,k)=0.
    151         enddo
    152       enddo
    153       if (prt_level.ge.1) print*,'7 OK convect8'
    154       do k=1,klev+1
    155          do ig=1,ngrid
    156             zw2(ig,k)=0.
    157             w_est(ig,k)=0.
    158             f_star(ig,k)=0.
    159             wa_moy(ig,k)=0.
     175      do l=1,klev
     176         do ig=1,ngrid
     177            if (alim_star_tot(ig) > 1.e-10 ) then
     178               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
     179            endif
    160180         enddo
    161181      enddo
    162 
    163       if (prt_level.ge.1) print*,'8 OK convect8'
    164       do ig=1,ngrid
    165          linter(ig)=1.
    166          lmix(ig)=1
    167          lmix_bis(ig)=2
    168          wmaxa(ig)=0.
    169       enddo
    170 
    171 !-----------------------------------------------------------------------------------
    172 !boucle de calcul de la vitesse verticale dans le thermique
    173 !-----------------------------------------------------------------------------------
    174       do l=1,klev-1
    175          do ig=1,ngrid
    176 
    177 
    178 
    179 ! Calcul dans la premiere couche active du thermique (ce qu'on teste
    180 ! en disant que la couche est instable et que w2 en bas de la couche
    181 ! est nulle.
    182 
    183             if (ztv(ig,l).gt.ztv(ig,l+1)  &
    184      &         .and.alim_star(ig,l).gt.1.e-10  &
    185      &         .and.zw2(ig,l).lt.1e-10) then
    186 
    187 
     182      alim_star_tot(:)=1.
     183
     184
     185!------------------------------------------------------------------------------
     186! Calcul dans la premiere couche
     187! On decide dans cette version que le thermique n'est actif que si la premiere
     188! couche est instable.
     189! Pourrait etre change si on veut que le thermiques puisse se déclencher
     190! dans une couche l>1
     191!------------------------------------------------------------------------------
     192do ig=1,ngrid
    188193! Le panache va prendre au debut les caracteristiques de l'air contenu
    189194! dans cette couche.
    190                ztla(ig,l)=zthl(ig,l)
    191                zqta(ig,l)=po(ig,l)
    192                zqla(ig,l)=zl(ig,l)
    193                f_star(ig,l+1)=alim_star(ig,l)
    194 
    195                zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
    196      &                     *(zlev(ig,l+1)-zlev(ig,l))  &
    197      &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
    198                w_est(ig,l+1)=zw2(ig,l+1)
     195    if (active(ig)) then
     196    ztla(ig,1)=zthl(ig,1)
     197    zqta(ig,1)=po(ig,1)
     198    zqla(ig,1)=zl(ig,1)
     199!cr: attention, prise en compte de f*(1)=1
     200    f_star(ig,2)=alim_star(ig,1)
     201    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
     202&                     *(zlev(ig,2)-zlev(ig,1))  &
     203&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
     204    w_est(ig,2)=zw2(ig,2)
     205    endif
     206enddo
    199207!
    200208
    201 
    202             else if ((zw2(ig,l).ge.1e-10).and.  &
    203      &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
    204 !estimation du detrainement a partir de la geometrie du pas precedent
    205 !tests sur la definition du detr
    206 !calcul de detr_star et entr_star
    207 
    208 
    209 
    210 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    211 ! FH le test miraculeux de Catherine ? Le bout du tunel ?
    212 !               w_est(ig,3)=zw2(ig,2)*  &
    213 !    &                   ((f_star(ig,2))**2)  &
    214 !    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
    215 !    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    216 !    &                   *(zlev(ig,3)-zlev(ig,2))
    217 !     if (l.gt.2) then
    218 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     209!==============================================================================
     210!boucle de calcul de la vitesse verticale dans le thermique
     211!==============================================================================
     212do l=2,klev-1
     213!==============================================================================
     214
     215
     216! On decide si le thermique est encore actif ou non
     217! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     218    do ig=1,ngrid
     219       active(ig)=active(ig) &
     220&                 .and. zw2(ig,l)>1.e-10 &
     221&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     222    enddo
    219223
    220224
    … …  
    222226! Premier calcul de la vitesse verticale a partir de la temperature
    223227! potentielle virtuelle
    224 
    225 ! FH CESTQUOI CA ????
    226 #define int1d2
    227 !#undef int1d2
    228 #ifdef int1d2
    229       if (l.ge.2) then
    230 #else
    231       if (l.gt.2) then
    232 #endif
    233 
    234       if (1.eq.1) then
    235           w_est(ig,3)=zw2(ig,2)* &
    236      &      ((f_star(ig,2))**2) &
    237      &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
    238      &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
    239 !     &      *1./3. &
    240      &      *(zlev(ig,3)-zlev(ig,2))
    241        endif
     228!     if (1.eq.1) then
     229!         w_est(ig,3)=zw2(ig,2)* &
     230!    &      ((f_star(ig,2))**2) &
     231!    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
     232!    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
     233!    &      *(zlev(ig,3)-zlev(ig,2))
     234!      endif
    242235
    243236
    244237!---------------------------------------------------------------------------
    245 !calcul de l entrainement et du detrainement lateral
     238! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     239! sans tenir compte du detrainement et de l'entrainement dans cette
     240! couche
     241! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     242! avant) a l'alimentation pour avoir un calcul plus propre
    246243!---------------------------------------------------------------------------
    247 !
    248 !test:estimation de ztva_new_est sans entrainement
    249 
    250                Tbef=ztla(ig,l-1)*zpspsk(ig,l)
    251                zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    252                qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
    253                qsatbef=MIN(0.5,qsatbef)
    254                zcor=1./(1.-retv*qsatbef)
    255                qsatbef=qsatbef*zcor
    256                Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
    257                if (Zsat) then
    258                qlbef=max(0.,zqta(ig,l-1)-qsatbef)
    259                DT = 0.5*RLvCp*qlbef
    260                do while (abs(DT).gt.DDT0)
    261                  Tbef=Tbef+DT
    262                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    263                  qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
    264                  qsatbef=MIN(0.5,qsatbef)
    265                  zcor=1./(1.-retv*qsatbef)
    266                  qsatbef=qsatbef*zcor
    267                  qlbef=zqta(ig,l-1)-qsatbef
    268 
    269                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    270                  zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
    271                  zcor=1./(1.-retv*qsatbef)
    272                  dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
    273                  num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
    274                  denom=1.+RLvCp*dqsat_dT
    275                  DT=num/denom
    276                enddo
    277                  zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
    278                endif
     244
     245     call thermcell_condens(ngrid,active, &
     246&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
     247
     248    do ig=1,ngrid
     249        if(active(ig)) then
    279250        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
     251        zta_est(ig,l)=ztva_est(ig,l)
    280252        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    281         zta_est(ig,l)=ztva_est(ig,l)
    282253        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    283254     &      -zqla_est(ig,l))-zqla_est(ig,l))
    … …  
    287258     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
    288259     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    289 !     &                   *1./3. &
    290260     &                   *(zlev(ig,l+1)-zlev(ig,l))
    291261             if (w_est(ig,l+1).lt.0.) then
    292262                w_est(ig,l+1)=zw2(ig,l)
    293263             endif
    294 !
    295 !calcul du detrainement
    296 !=======================
    297 
    298 !CR:on vire les modifs
    299          if (iflag_thermals_ed==0) then
    300 
    301 ! Modifications du calcul du detrainement.
    302 ! Dans la version de la these de Catherine, on passe brusquement
    303 ! de la version seche a la version nuageuse pour le detrainement
    304 ! ce qui peut occasioner des oscillations.
    305 ! dans la nouvelle version, on commence par calculer un detrainement sec.
    306 ! Puis un autre en cas de nuages.
    307 ! Puis on combine les deux lineairement en fonction de la quantite d'eau.
    308 
    309 #define int1d3
    310 !#undef int1d3
    311 #define RIO_TH
    312 #ifdef RIO_TH
    313 !1. Cas non nuageux
    314 ! 1.1 on est sous le zmax_sec et w croit
    315           if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
    316      &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
    317 #ifdef int1d3
    318      &       (zqla_est(ig,l).lt.1.e-10)) then
    319 #else
    320      &       (zqla(ig,l-1).lt.1.e-10)) then
    321 #endif
    322              detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
    323      &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
    324      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
    325      &       /(r_aspect*zmax_sec(ig)))
    326              detr_stara(ig,l)=detr_star(ig,l)
    327 
    328        if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
    329 
    330 ! 1.2 on est sous le zmax_sec et w decroit
    331           else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
    332 #ifdef int1d3
    333      &            (zqla_est(ig,l).lt.1.e-10)) then
    334 #else
    335      &            (zqla(ig,l-1).lt.1.e-10)) then
    336 #endif
    337              detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
    338      &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
    339      &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
    340      &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
    341      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
    342      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
    343      &       *((zmax_sec(ig)-zlev(ig,l))/  &
    344      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
    345              detr_starb(ig,l)=detr_star(ig,l)
    346 
    347         if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
    348 
     264       endif
     265    enddo
     266
     267!-------------------------------------------------
     268!calcul des taux d'entrainement et de detrainement
     269!-------------------------------------------------
     270
     271     do ig=1,ngrid
     272        if (active(ig)) then
     273          zdz=zlev(ig,l+1)-zlev(ig,l)
     274          zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     275          zfact=fact_gamma/(1.+fact_gamma)
     276 
     277! estimation de la fraction couverte par les thermiques
     278          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
     279
     280!calcul de la soumission papier
     281          if (1.eq.1) then
     282             fact_epsilon=0.0007
     283!             zfact=0.9/(1.+0.9)
     284             zfact=0.3
     285             fact_gamma=0.7
     286             fact_gamma2=0.6
     287             expa=0.25
     288! Calcul  du taux d'entrainement entr_star (epsilon)
     289           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
     290     &     zbuoy/w_est(ig,l+1) )&
     291!- fact_epsilon/zalpha**0.25  ) &
     292     &     +0.000 )
     293
     294!           entr_star(ig,l)=f_star(ig,l)*zdz * (  1./3 * MAX(0.,  &     
     295!     &     zbuoy/w_est(ig,l+1) - 1./zalpha**0.25  ) &
     296!     &     +0.000 )
     297! Calcul du taux de detrainement detr_star (delta)
     298!           if (zqla_est(ig,l).lt.1.e-10) then
     299!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     300!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
     301!     &     +0.0006 )
     302!           else
     303!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     304!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
     305!     &     +0.002 )
     306!           endif
     307           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     308     &     fact_gamma2 * MAX(0., &
     309!fact_epsilon/zalpha**0.25
     310     &      -zbuoy/w_est(ig,l+1) )       &
     311!     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
     312!test
     313     &     +  0.006*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &   
     314     &     +0.0000 )
    349315          else
    350316
    351 ! 1.3 dans les autres cas
    352              detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
    353      &                      *(zlev(ig,l+1)-zlev(ig,l))
    354              detr_starc(ig,l)=detr_star(ig,l)
    355 
    356         if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
    357              
     317! Calcul  du taux d'entrainement entr_star (epsilon)
     318           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
     319     &     zbuoy/w_est(ig,l+1) - fact_epsilon  ) &
     320     &     +0.0000 )
     321
     322! Calcul du taux de detrainement detr_star (delta)
     323           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     324     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
     325     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
     326     &     +0.0000 )
     327
    358328          endif
    359 
    360 #else
    361 
    362 ! 1.1 on est sous le zmax_sec et w croit
    363           if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
    364      &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
    365              detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
    366      &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
    367      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
    368      &       /(r_aspect*zmax_sec(ig)))
    369 
    370        if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
    371 
    372 ! 1.2 on est sous le zmax_sec et w decroit
    373           else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
    374              detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
    375      &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
    376      &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
    377      &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
    378      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
    379      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
    380      &       *((zmax_sec(ig)-zlev(ig,l))/  &
    381      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
    382        if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
    383 
    384           else
    385              detr_star=0.
    386           endif
    387 
    388 ! 1.3 dans les autres cas
    389           detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
    390      &                      *(zlev(ig,l+1)-zlev(ig,l))
    391 
    392           coefc=min(zqla(ig,l-1)/1.e-3,1.)
    393           if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
    394           coefc=1.
    395 ! il semble qu'il soit important de baser le calcul sur
    396 ! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
    397           detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
    398 
    399         if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
    400 
    401 #endif
    402 
    403 
    404         if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
    405 !IM 730508 beg
    406 !        if(itap.GE.7200) THEN
    407 !         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
    408 !        endif
    409 !IM 730508 end
    410          
    411          zqla0(ig,l)=zqla_est(ig,l)
    412          detr_star0(ig,l)=detr_star(ig,l)
    413 !IM 060508 beg
    414 !         if(detr_star(ig,l).GT.1.) THEN
    415 !          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
    416 !   &      ,detr_starc(ig,l),coefc
    417 !         endif
    418 !IM 060508 end
    419 !IM 160508 beg
    420 !IM 160508       IF (f0(ig).NE.0.) THEN
    421            detr_star(ig,l)=detr_star(ig,l)/f0(ig)
    422 !IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
    423 !IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
    424 !IM 160508       ELSE
    425 !IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
    426 !IM 160508       ENDIF
    427 !IM 160508 end
    428 !IM 060508 beg
    429 !        if(detr_star(ig,l).GT.1.) THEN
    430 !         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
    431 !   &     float(1)/f0(ig)
    432 !        endif
    433 !IM 060508 end
    434         if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
    435 !
    436 !calcul de entr_star
    437 
    438 ! #undef test2
    439 ! #ifdef test2
    440 ! La version test2 destabilise beaucoup le modele.
    441 ! Il semble donc que ca aide d'avoir un entrainement important sous
    442 ! le nuage.
    443 !         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
    444 !          entr_star(ig,l)=0.4*detr_star(ig,l)
    445 !         else
    446 !          entr_star(ig,l)=0.
    447 !         endif
    448 ! #else
    449 !
    450 ! Deplacement du calcul de entr_star pour eviter d'avoir aussi
    451 ! entr_star > fstar.
    452 ! Redeplacer suite a la transformation du cas detr>f
    453 ! FH
    454 
    455         if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
    456 #define int1d
    457 !FH 070508 #define int1d4
    458 !#undef int1d4
    459 ! L'option int1d4 correspond au choix dans le cas ou le detrainement
    460 ! devient trop grand.
    461 
    462 #ifdef int1d
    463 
    464 #ifdef int1d4
    465 #else
    466        detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
    467 !FH 070508 plus
    468        detr_star(ig,l)=min(detr_star(ig,l),1.)
    469 #endif
    470 
    471        entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
    472 
    473         if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
    474 #ifdef int1d4
    475 ! Si le detrainement excede le flux en bas + l'entrainement, le thermique
    476 ! doit disparaitre.
    477        if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
    478           detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
    479           f_star(ig,l+1)=0.
    480           linter(ig)=l+1
    481           zw2(ig,l+1)=-1.e-10
    482        endif
    483 #endif
    484 
    485 
    486 #else
    487 
    488         if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
    489         if(l.gt.lalim(ig)) then
    490          entr_star(ig,l)=0.4*detr_star(ig,l)
    491         else
    492 
    493 ! FH :
    494 ! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
    495 ! en haut de la couche d'alimentation.
    496 ! A remettre en questoin a la premiere occasion mais ca peut aider a
    497 ! ecrire un code robuste.
    498 ! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
    499 ! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
    500 ! d'alimentation, ce qui n'est pas forcement heureux.
    501 
    502         if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
    503 #undef pre_int1c
    504 #ifdef pre_int1c
    505          entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
    506          detr_star(ig,l)=entr_star(ig,l)
    507 #else
    508          entr_star(ig,l)=0.
    509 #endif
    510 
     329!endif choix du calcul de E* et D*
     330
     331!cr test
     332!           entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l)     
     333
     334! Prise en compte de la fraction
     335!      detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5))
     336
     337! En dessous de lalim, on prend le max de alim_star et entr_star pour
     338! alim_star et 0 sinon
     339        if (l.lt.lalim(ig)) then
     340          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
     341          entr_star(ig,l)=0.
    511342        endif
    512343
    513 #endif
    514 
    515         if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
    516         entr_star1(ig,l)=entr_star(ig,l)
    517         detr_star1(ig,l)=detr_star(ig,l)
    518 !
    519 
    520 #ifdef int1d
    521 #else
    522         if (detr_star(ig,l).gt.f_star(ig,l)) then
    523 
    524 !  Ce test est là entre autres parce qu'on passe par des valeurs
    525 !  delirantes de detr_star.
    526 !  ca vaut sans doute le coup de verifier pourquoi.
    527 
    528            detr_star(ig,l)=f_star(ig,l)
    529 #ifdef pre_int1c
    530            if (l.gt.lalim(ig)+1) then
    531                entr_star(ig,l)=0.
    532                alim_star(ig,l)=0.
    533 ! FH ajout pour forcer a stoper le thermique juste sous le sommet
    534 ! de la couche (voir calcul de finter)
    535                zw2(ig,l+1)=-1.e-10
    536                linter(ig)=l+1
    537             else
    538                entr_star(ig,l)=0.4*detr_star(ig,l)
    539             endif
    540 #else
    541            entr_star(ig,l)=0.4*detr_star(ig,l)
    542 #endif
    543         endif
    544 #endif
    545 
    546       else !l > 2
    547          detr_star(ig,l)=0.
    548          entr_star(ig,l)=0.
    549       endif
    550 
    551         entr_star2(ig,l)=entr_star(ig,l)
    552         detr_star2(ig,l)=detr_star(ig,l)
    553         if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
    554 
    555        endif  ! iflag_thermals_ed==0
    556 
    557 !CR:nvlle def de entr_star et detr_star
    558       if (iflag_thermals_ed>=1) then
    559 !      if (l.lt.lalim(ig)) then
    560 !      if (l.lt.2) then
    561 !        entr_star(ig,l)=0.
    562 !        detr_star(ig,l)=0.
    563 !      else
    564 !      if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then
    565 !         entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    566 !      else
    567 !         entr_star(ig,l)=  &
    568 !     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
    569 !     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
    570 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    571 
    572  
    573          entr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &         
    574      &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
    575      &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)   &
    576      &                *(zlev(ig,l+1)-zlev(ig,l))) &
    577      &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    578 
    579         if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
    580             alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
    581             lalim(ig)=lmix_bis(ig)
    582             if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
    583         endif
    584 
    585         if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
    586 !        c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
    587          c2(ig,l)=0.001
    588          detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
    589      &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
    590      &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
    591      &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
    592      &                *(zlev(ig,l+1)-zlev(ig,l)))                    &
    593      &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    594 
    595        else
    596 !         c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
    597           c2(ig,l)=0.003
    598 
    599          detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
    600      &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
    601      &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
    602      &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
    603      &                *(zlev(ig,l+1)-zlev(ig,l))) &
    604      &                +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    605        endif
    606          
    607            
    608 !        detr_star(ig,l)=detr_star(ig,l)*3.
    609 !        if (l.lt.lalim(ig)) then
    610 !          entr_star(ig,l)=0.
    611 !        endif
    612 !        if (l.lt.2) then
    613 !          entr_star(ig,l)=0.
    614 !          detr_star(ig,l)=0.
    615 !        endif
    616 
    617 
    618 !      endif
    619 !      else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
    620 !      entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1))        &
    621 !     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
    622 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    623 !      detr_star(ig,l)=0.002*f_star(ig,l)                         &
    624 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    625 !      else
    626 !      entr_star(ig,l)=0.001*f_star(ig,l)                         &
    627 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    628 !      detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1))       &
    629 !     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))       &
    630 !     &                *(zlev(ig,l+1)-zlev(ig,l))                      &
    631 !     &                +0.002*f_star(ig,l)                             &
    632 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    633 !      endif
    634 
    635       endif   ! iflag_thermals_ed==1
    636 
    637 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    638 ! FH inutile si on conserve comme on l'a fait plus haut entr=detr
    639 ! dans la couche d'alimentation
    640 !pas d entrainement dans la couche alim
    641 !      if ((l.le.lalim(ig))) then
    642 !           entr_star(ig,l)=0.
    643 !      endif
    644 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    645 !
    646 !prise en compte du detrainement et de l entrainement dans le calcul du flux
    647 
     344! Calcul du flux montant normalise
    648345      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    649346     &              -detr_star(ig,l)
    650347
    651 !test sur le signe de f_star
    652         if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
    653        if (f_star(ig,l+1).gt.1.e-10) then
     348      endif
     349   enddo
     350
    654351!----------------------------------------------------------------------------
    655352!calcul de la vitesse verticale en melangeant Tl et qt du thermique
    656353!---------------------------------------------------------------------------
    657 !
    658        Zsat=.false.
    659        ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
     354   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
     355   do ig=1,ngrid
     356       if (activetmp(ig)) then
     357           Zsat=.false.
     358           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    660359     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
    661360     &            /(f_star(ig,l+1)+detr_star(ig,l))
    662 !
    663        zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
     361           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
    664362     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
    665363     &            /(f_star(ig,l+1)+detr_star(ig,l))
    666 ! 
    667                Tbef=ztla(ig,l)*zpspsk(ig,l)
    668                zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    669                qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
    670                qsatbef=MIN(0.5,qsatbef)
    671                zcor=1./(1.-retv*qsatbef)
    672                qsatbef=qsatbef*zcor
    673                Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
    674                if (Zsat) then
    675                qlbef=max(0.,zqta(ig,l)-qsatbef)
    676                DT = 0.5*RLvCp*qlbef
    677                do while (abs(DT).gt.DDT0)
    678                  Tbef=Tbef+DT
    679                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    680                  qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
    681                  qsatbef=MIN(0.5,qsatbef)
    682                  zcor=1./(1.-retv*qsatbef)
    683                  qsatbef=qsatbef*zcor
    684                  qlbef=zqta(ig,l)-qsatbef
    685 
    686                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    687                  zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
    688                  zcor=1./(1.-retv*qsatbef)
    689                  dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
    690                  num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
    691                  denom=1.+RLvCp*dqsat_dT
    692                  DT=num/denom
    693               enddo
    694                  zqla(ig,l) = max(0.,qlbef)
    695               endif
    696 !   
     364
     365        endif
     366    enddo
     367
     368   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
     369
     370
     371   do ig=1,ngrid
     372      if (activetmp(ig)) then
    697373        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
    698374! on ecrit de maniere conservative (sat ou non)
    … …  
    707383!on ecrit zqsat
    708384           zqsatth(ig,l)=qsatbef 
    709 !calcul de vitesse
    710            zw2(ig,l+1)=zw2(ig,l)*  &
    711      &                 ((f_star(ig,l))**2)  &
    712 !  Tests de Catherine
    713 !     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
    714      &      /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ &
    715      &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    716      &                 *fact_gamma &
    717      &                 *(zlev(ig,l+1)-zlev(ig,l))
    718 !prise en compte des forces de pression que qd flottabilité<0
    719 !              zw2(ig,l+1)=zw2(ig,l)*  &
    720 !     &            1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &       
    721 !     &                 (f_star(ig,l))**2 &
    722 !     &                 /(f_star(ig,l)+entr_star(ig,l))**2+ &
    723 !     &                 (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+  &       
    724 !     &                 /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ &
    725 !     &                 /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ &
    726 !     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    727 !     &                 *1./3. &
     385
     386!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     387!          zw2(ig,l+1)=&
     388!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
     389!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
     390!     &                 *1./(1.+fact_gamma) &
    728391!     &                 *(zlev(ig,l+1)-zlev(ig,l))
    729          
    730 !        write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
    731 !     &              -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), &
    732 !     &               2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    733 
    734  
    735 !             zw2(ig,l+1)=zw2(ig,l)*  &
    736 !     &                 (2.-2.*entr_star(ig,l)/f_star(ig,l)) & 
    737 !     &                 -zw2(ig,l-1)+  &       
    738 !     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    739 !     &                 *1./3. &
    740 !     &                 *(zlev(ig,l+1)-zlev(ig,l))             
    741 
    742             endif
    743         endif
     392!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     393! La meme en plus modulaire :
     394           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     395           zdz=zlev(ig,l+1)-zlev(ig,l)
     396
     397
     398           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
     399
     400if (1==0) then
     401           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
     402           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
     403           zw2(ig,l+1)=zw2modif+zdw2
     404else
     405! Tentative de reecriture de l'equation de w2. A reprendre ...
     406!           zdw2=2*zdz*zbuoy
     407!           zw2modif=zw2(ig,l)*(1.-2.*zdz*(zeps+fact_epsilon))
     408!!!!!      zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l))
     409!formulation Arnaud
     410!           zw2fact=zbuoy*zalpha**expa/fact_epsilon
     411!           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa*(1+fact_gamma))*zdz) &
     412!      &    +zw2fact
     413           if (zbuoy.gt.1.e-10) then
     414           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma-zfact)
     415           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
     416      &    +zw2fact
     417           else
     418           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma)
     419           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
     420      &    +zw2fact
     421
     422           endif
     423
     424endif
     425!           zw2(ig,l+1)=zw2modif+zdw2
     426      endif
     427   enddo
     428
    744429        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
    745430!
     431!---------------------------------------------------------------------------
    746432!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    747 
     433!---------------------------------------------------------------------------
     434
     435   do ig=1,ngrid
    748436            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    749437!               stop'On tombe sur le cas particulier de thermcell_dry'
    … …  
    753441            endif
    754442
    755 !        if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
    756443        if (zw2(ig,l+1).lt.0.) then
    757444           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    … …  
    771458            wmaxa(ig)=wa_moy(ig,l+1)
    772459        endif
    773         enddo
     460   enddo
     461
     462!=========================================================================
     463! FIN DE LA BOUCLE VERTICALE
    774464      enddo
    775 
    776 !on remplace a* par e* ds premiere couche
    777 !      if (iflag_thermals_ed.ge.1) then
    778 !       do ig=1,ngrid
    779 !       do l=2,klev
    780 !          if (l.lt.lalim(ig)) then
    781 !             alim_star(ig,l)=entr_star(ig,l)
    782 !          endif
    783 !       enddo
    784 !       enddo
    785 !       do ig=1,ngrid
    786 !          lalim(ig)=lmix_bis(ig)
    787 !       enddo
    788 !      endif
    789        if (iflag_thermals_ed.ge.1) then
    790           do ig=1,ngrid
    791              do l=2,lalim(ig)
    792                 alim_star(ig,l)=entr_star(ig,l)
    793                 entr_star(ig,l)=0.
    794              enddo
    795            enddo
    796        endif
     465!=========================================================================
     466
     467!on recalcule alim_star_tot
     468       do ig=1,ngrid
     469          alim_star_tot(ig)=0.
     470       enddo
     471       do ig=1,ngrid
     472          do l=1,lalim(ig)-1
     473          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     474          enddo
     475       enddo
     476       
     477
    797478        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    798479
Note: See TracChangeset for help on using the changeset viewer.