Ignore:
Timestamp:
Jul 1, 2010, 11:02:53 AM (14 years ago)
Author:
Laurent Fairhead
Message:

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

Location:
LMDZ4/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk

  • LMDZ4/trunk/libf/phylmd/thermcell_plume.F90

    r1057 r1403  
     1!
     2! $Id$
     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,zdrag
    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)
    90 
     87      REAL a1,m
     88
     89      REAL zw2fact,expa
    9190      Zsat=.false.
    9291! Initialisation
    9392      RLvCp = RLVTT/RCPD
    9493     
    95       if (iflag_thermals_ed==0) then
    96          fact_gamma=1.
    97          fact_epsilon=1.
    98       else if (iflag_thermals_ed==1)  then
    99          fact_gamma=1.
    100          fact_epsilon=1.
    101       else if (iflag_thermals_ed==2)  then
    102          fact_gamma=1.
    103          fact_epsilon=2.
    104       endif
    105 
    106       do l=1,klev
     94     
     95         fact_epsilon=0.002
     96         a1=2./3.
     97         fact_gamma=0.9
     98         zfact=fact_gamma/(1+fact_gamma)
     99         fact_gamma2=zfact
     100         expa=0.
     101     
     102
     103! Initialisations des variables reeles
     104if (1==1) then
     105      ztva(:,:)=ztv(:,:)
     106      ztva_est(:,:)=ztva(:,:)
     107      ztla(:,:)=zthl(:,:)
     108      zqta(:,:)=po(:,:)
     109      zha(:,:) = ztva(:,:)
     110else
     111      ztva(:,:)=0.
     112      ztva_est(:,:)=0.
     113      ztla(:,:)=0.
     114      zqta(:,:)=0.
     115      zha(:,:) =0.
     116endif
     117
     118      zqla_est(:,:)=0.
     119      zqsatth(:,:)=0.
     120      zqla(:,:)=0.
     121      detr_star(:,:)=0.
     122      entr_star(:,:)=0.
     123      alim_star(:,:)=0.
     124      alim_star_tot(:)=0.
     125      csc(:,:)=0.
     126      detr(:,:)=0.
     127      entr(:,:)=0.
     128      zw2(:,:)=0.
     129      w_est(:,:)=0.
     130      f_star(:,:)=0.
     131      wa_moy(:,:)=0.
     132      linter(:)=1.
     133      linter(:)=1.
     134
     135! Initialisation des variables entieres
     136      lmix(:)=1
     137      lmix_bis(:)=2
     138      wmaxa(:)=0.
     139      lalim(:)=1
     140
     141!-------------------------------------------------------------------------
     142! On ne considere comme actif que les colonnes dont les deux premieres
     143! couches sont instables.
     144!-------------------------------------------------------------------------
     145      active(:)=ztv(:,1)>ztv(:,2)
     146
     147!-------------------------------------------------------------------------
     148! Definition de l'alimentation a l'origine dans thermcell_init
     149!-------------------------------------------------------------------------
     150      do l=1,klev-1
    107151         do ig=1,ngrid
    108             zqla_est(ig,l)=0.
    109             ztva_est(ig,l)=ztva(ig,l)
    110             zqsatth(ig,l)=0.
     152            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
     153               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
     154     &                       *sqrt(zlev(ig,l+1))
     155               lalim(:)=l+1
     156               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     157            endif
    111158         enddo
    112159      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.
     160      do l=1,klev
     161         do ig=1,ngrid
     162            if (alim_star_tot(ig) > 1.e-10 ) then
     163               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
     164            endif
    160165         enddo
    161166      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 
     167      alim_star_tot(:)=1.
     168
     169
     170!------------------------------------------------------------------------------
     171! Calcul dans la premiere couche
     172! On decide dans cette version que le thermique n'est actif que si la premiere
     173! couche est instable.
     174! Pourrait etre change si on veut que le thermiques puisse se déclencher
     175! dans une couche l>1
     176!------------------------------------------------------------------------------
     177do ig=1,ngrid
    188178! Le panache va prendre au debut les caracteristiques de l'air contenu
    189179! 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)
     180    if (active(ig)) then
     181    ztla(ig,1)=zthl(ig,1)
     182    zqta(ig,1)=po(ig,1)
     183    zqla(ig,1)=zl(ig,1)
     184!cr: attention, prise en compte de f*(1)=1
     185    f_star(ig,2)=alim_star(ig,1)
     186    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
     187&                     *(zlev(ig,2)-zlev(ig,1))  &
     188&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
     189    w_est(ig,2)=zw2(ig,2)
     190    endif
     191enddo
    199192!
    200193
    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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     194!==============================================================================
     195!boucle de calcul de la vitesse verticale dans le thermique
     196!==============================================================================
     197do l=2,klev-1
     198!==============================================================================
     199
     200
     201! On decide si le thermique est encore actif ou non
     202! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     203    do ig=1,ngrid
     204       active(ig)=active(ig) &
     205&                 .and. zw2(ig,l)>1.e-10 &
     206&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     207    enddo
    219208
    220209
    … …  
    222211! Premier calcul de la vitesse verticale a partir de la temperature
    223212! 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
    242 
    243 
    244 !---------------------------------------------------------------------------
    245 !calcul de l entrainement et du detrainement lateral
    246 !---------------------------------------------------------------------------
    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
     213!     if (1.eq.1) then
     214!         w_est(ig,3)=zw2(ig,2)* &
     215!    &      ((f_star(ig,2))**2) &
     216!    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
     217!    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
     218!    &      *(zlev(ig,3)-zlev(ig,2))
     219!      endif
     220
     221
     222!---------------------------------------------------------------------------
     223! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     224! sans tenir compte du detrainement et de l'entrainement dans cette
     225! couche
     226! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     227! avant) a l'alimentation pour avoir un calcul plus propre
     228!---------------------------------------------------------------------------
     229
     230     call thermcell_condens(ngrid,active, &
     231&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
     232
     233    do ig=1,ngrid
     234        if(active(ig)) then
    279235        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
     236        zta_est(ig,l)=ztva_est(ig,l)
    280237        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    281         zta_est(ig,l)=ztva_est(ig,l)
    282238        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    283239     &      -zqla_est(ig,l))-zqla_est(ig,l))
    284240
     241         if (1.eq.0) then
     242!calcul de w_est sans prendre en compte le drag
    285243             w_est(ig,l+1)=zw2(ig,l)*  &
    286244     &                   ((f_star(ig,l))**2)  &
    287245     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
    288246     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    289 !     &                   *1./3. &
    290247     &                   *(zlev(ig,l+1)-zlev(ig,l))
     248         else
     249
     250           zdz=zlev(ig,l+1)-zlev(ig,l)
     251           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
     252           zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     253           zdrag=fact_epsilon/(zalpha**expa)
     254           zw2fact=zbuoy/zdrag*a1
     255           w_est(ig,l+1)=(w_est(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
     256      &    +zw2fact
     257
     258         endif
     259
    291260             if (w_est(ig,l+1).lt.0.) then
    292261                w_est(ig,l+1)=zw2(ig,l)
    293262             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 
    349           else
    350 
    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              
    358           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)
     263       endif
     264    enddo
     265
     266!-------------------------------------------------
     267!calcul des taux d'entrainement et de detrainement
     268!-------------------------------------------------
     269
     270     do ig=1,ngrid
     271        if (active(ig)) then
     272          zdz=zlev(ig,l+1)-zlev(ig,l)
     273          zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     274 
     275! estimation de la fraction couverte par les thermiques
     276          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
     277
     278!calcul de la soumission papier
     279! Calcul  du taux d'entrainement entr_star (epsilon)
     280           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
     281     &     a1*zbuoy/w_est(ig,l+1) &
     282     &     - fact_epsilon/zalpha**expa  ) &
     283     &     +0. )
     284           
     285!calcul du taux de detrainment (delta)
     286!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     287!     &      MAX(1.e-3, &
     288!     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
     289!     &      +0.01*(max(zqta(ig,l-1)-po(ig,l),0.)/(po(ig,l))/(w_est(ig,l+1)))**0.5    &   
     290!     &     +0. ))
     291
     292          m=0.5
     293
     294          detr_star(ig,l)=1.*f_star(ig,l)*zdz *                    &
     295    &     MAX(5.e-4,-fact_gamma2*a1*(1./w_est(ig,l+1))*((1.-(1.-m)/(1.+70*zqta(ig,l-1)))*zbuoy        &
     296    &     -40*(1.-m)*(max(zqta(ig,l-1)-po(ig,l),0.))/(1.+70*zqta(ig,l-1)) )   )
     297
     298!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     299!     &      MAX(0.0, &
     300!     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
     301!     &      +20*(max(zqta(ig,l-1)-po(ig,l),0.))**1*(zalpha/w_est(ig,l+1))**0.5    &   
     302!     &     +0. ))
     303
     304
     305! En dessous de lalim, on prend le max de alim_star et entr_star pour
     306! alim_star et 0 sinon
     307        if (l.lt.lalim(ig)) then
     308          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
     309          entr_star(ig,l)=0.
     310        endif
     311
     312!attention test
     313!        if (detr_star(ig,l).gt.(f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l))) then       
     314!            detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
    408315!        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 
    511         endif
    512 
    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 
     316! Calcul du flux montant normalise
    648317      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    649318     &              -detr_star(ig,l)
    650319
    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
     320      endif
     321   enddo
     322
    654323!----------------------------------------------------------------------------
    655324!calcul de la vitesse verticale en melangeant Tl et qt du thermique
    656325!---------------------------------------------------------------------------
    657 !
    658        Zsat=.false.
    659        ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
     326   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
     327   do ig=1,ngrid
     328       if (activetmp(ig)) then
     329           Zsat=.false.
     330           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    660331     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
    661332     &            /(f_star(ig,l+1)+detr_star(ig,l))
    662 !
    663        zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
     333           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
    664334     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
    665335     &            /(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 !   
     336
     337        endif
     338    enddo
     339
     340   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
     341
     342
     343   do ig=1,ngrid
     344      if (activetmp(ig)) then
    697345        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
    698346! on ecrit de maniere conservative (sat ou non)
    … …  
    707355!on ecrit zqsat
    708356           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. &
     357
     358!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     359!          zw2(ig,l+1)=&
     360!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
     361!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
     362!     &                 *1./(1.+fact_gamma) &
    728363!     &                 *(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
     364!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     365! La meme en plus modulaire :
     366           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     367           zdz=zlev(ig,l+1)-zlev(ig,l)
     368
     369
     370           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
     371
     372if (1==0) then
     373           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
     374           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
     375           zw2(ig,l+1)=zw2modif+zdw2
     376else
     377           zdrag=fact_epsilon/(zalpha**expa)
     378           zw2fact=zbuoy/zdrag*a1
     379           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
     380      &    +zw2fact
     381
     382
     383endif
     384
     385      endif
     386   enddo
     387
    744388        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
    745389!
     390!---------------------------------------------------------------------------
    746391!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    747 
     392!---------------------------------------------------------------------------
     393
     394   do ig=1,ngrid
    748395            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    749396!               stop'On tombe sur le cas particulier de thermcell_dry'
    … …  
    753400            endif
    754401
    755 !        if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
    756402        if (zw2(ig,l+1).lt.0.) then
    757403           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    … …  
    771417            wmaxa(ig)=wa_moy(ig,l+1)
    772418        endif
    773         enddo
     419   enddo
     420
     421!=========================================================================
     422! FIN DE LA BOUCLE VERTICALE
    774423      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
     424!=========================================================================
     425
     426!on recalcule alim_star_tot
     427       do ig=1,ngrid
     428          alim_star_tot(ig)=0.
     429       enddo
     430       do ig=1,ngrid
     431          do l=1,lalim(ig)-1
     432          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     433          enddo
     434       enddo
     435       
     436
    797437        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    798438
    799 !     print*,'thermcell_plume OK'
    800439
    801440      return
    802441      end
     442
     443!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     444!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     445!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     446!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     447!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     448!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     449 SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
     450&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
     451&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     452&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     453&           ,lev_out,lunout1,igout)
     454
     455!--------------------------------------------------------------------------
     456!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
     457! Version conforme a l'article de Rio et al. 2010.
     458! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
     459!--------------------------------------------------------------------------
     460
     461      IMPLICIT NONE
     462
     463#include "YOMCST.h"
     464#include "YOETHF.h"
     465#include "FCTTRE.h"
     466#include "iniprint.h"
     467#include "thermcell.h"
     468
     469      INTEGER itap
     470      INTEGER lunout1,igout
     471      INTEGER ngrid,klev
     472      REAL ptimestep
     473      REAL ztv(ngrid,klev)
     474      REAL zthl(ngrid,klev)
     475      REAL po(ngrid,klev)
     476      REAL zl(ngrid,klev)
     477      REAL rhobarz(ngrid,klev)
     478      REAL zlev(ngrid,klev+1)
     479      REAL pplev(ngrid,klev+1)
     480      REAL pphi(ngrid,klev)
     481      REAL zpspsk(ngrid,klev)
     482      REAL alim_star(ngrid,klev)
     483      REAL f0(ngrid)
     484      INTEGER lalim(ngrid)
     485      integer lev_out                           ! niveau pour les print
     486   
     487      real alim_star_tot(ngrid)
     488
     489      REAL ztva(ngrid,klev)
     490      REAL ztla(ngrid,klev)
     491      REAL zqla(ngrid,klev)
     492      REAL zqta(ngrid,klev)
     493      REAL zha(ngrid,klev)
     494
     495      REAL detr_star(ngrid,klev)
     496      REAL coefc
     497      REAL entr_star(ngrid,klev)
     498      REAL detr(ngrid,klev)
     499      REAL entr(ngrid,klev)
     500
     501      REAL csc(ngrid,klev)
     502
     503      REAL zw2(ngrid,klev+1)
     504      REAL w_est(ngrid,klev+1)
     505      REAL f_star(ngrid,klev+1)
     506      REAL wa_moy(ngrid,klev+1)
     507
     508      REAL ztva_est(ngrid,klev)
     509      REAL zqla_est(ngrid,klev)
     510      REAL zqsatth(ngrid,klev)
     511      REAL zta_est(ngrid,klev)
     512      REAL ztemp(ngrid),zqsat(ngrid)
     513      REAL zdw2
     514      REAL zw2modif
     515      REAL zw2fact
     516      REAL zeps(ngrid,klev)
     517
     518      REAL linter(ngrid)
     519      INTEGER lmix(ngrid)
     520      INTEGER lmix_bis(ngrid)
     521      REAL    wmaxa(ngrid)
     522
     523      INTEGER ig,l,k
     524
     525      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
     526      real zbuoybis
     527      real zcor,zdelta,zcvm5,qlbef,zdz2
     528      real betalpha,zbetalpha
     529      real eps, afact
     530      REAL REPS,RLvCp,DDT0
     531      PARAMETER (DDT0=.01)
     532      logical Zsat
     533      LOGICAL active(ngrid),activetmp(ngrid)
     534      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
     535      REAL c2(ngrid,klev)
     536      Zsat=.false.
     537! Initialisation
     538
     539      RLvCp = RLVTT/RCPD
     540      fact_epsilon=0.002
     541      betalpha=0.9
     542      afact=2./3.           
     543
     544      zbetalpha=betalpha/(1.+betalpha)
     545
     546
     547! Initialisations des variables reeles
     548if (1==0) then
     549      ztva(:,:)=ztv(:,:)
     550      ztva_est(:,:)=ztva(:,:)
     551      ztla(:,:)=zthl(:,:)
     552      zqta(:,:)=po(:,:)
     553      zha(:,:) = ztva(:,:)
     554else
     555      ztva(:,:)=0.
     556      ztva_est(:,:)=0.
     557      ztla(:,:)=0.
     558      zqta(:,:)=0.
     559      zha(:,:) =0.
     560endif
     561
     562      zqla_est(:,:)=0.
     563      zqsatth(:,:)=0.
     564      zqla(:,:)=0.
     565      detr_star(:,:)=0.
     566      entr_star(:,:)=0.
     567      alim_star(:,:)=0.
     568      alim_star_tot(:)=0.
     569      csc(:,:)=0.
     570      detr(:,:)=0.
     571      entr(:,:)=0.
     572      zw2(:,:)=0.
     573      zbuoy(:,:)=0.
     574      gamma(:,:)=0.
     575      zeps(:,:)=0.
     576      w_est(:,:)=0.
     577      f_star(:,:)=0.
     578      wa_moy(:,:)=0.
     579      linter(:)=1.
     580!     linter(:)=1.
     581! Initialisation des variables entieres
     582      lmix(:)=1
     583      lmix_bis(:)=2
     584      wmaxa(:)=0.
     585      lalim(:)=1
     586
     587
     588!-------------------------------------------------------------------------
     589! On ne considere comme actif que les colonnes dont les deux premieres
     590! couches sont instables.
     591!-------------------------------------------------------------------------
     592      active(:)=ztv(:,1)>ztv(:,2)
     593
     594!-------------------------------------------------------------------------
     595! Definition de l'alimentation a l'origine dans thermcell_init
     596!-------------------------------------------------------------------------
     597      do l=1,klev-1
     598         do ig=1,ngrid
     599            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
     600               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
     601     &                       *sqrt(zlev(ig,l+1))
     602               lalim(ig)=l+1
     603               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     604            endif
     605         enddo
     606      enddo
     607      do l=1,klev
     608         do ig=1,ngrid
     609            if (alim_star_tot(ig) > 1.e-10 ) then
     610               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
     611            endif
     612         enddo
     613      enddo
     614      alim_star_tot(:)=1.
     615
     616
     617
     618!------------------------------------------------------------------------------
     619! Calcul dans la premiere couche
     620! On decide dans cette version que le thermique n'est actif que si la premiere
     621! couche est instable.
     622! Pourrait etre change si on veut que le thermiques puisse se déclencher
     623! dans une couche l>1
     624!------------------------------------------------------------------------------
     625do ig=1,ngrid
     626! Le panache va prendre au debut les caracteristiques de l'air contenu
     627! dans cette couche.
     628    if (active(ig)) then
     629    ztla(ig,1)=zthl(ig,1)
     630    zqta(ig,1)=po(ig,1)
     631    zqla(ig,1)=zl(ig,1)
     632!cr: attention, prise en compte de f*(1)=1
     633    f_star(ig,2)=alim_star(ig,1)
     634    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
     635&                     *(zlev(ig,2)-zlev(ig,1))  &
     636&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
     637    w_est(ig,2)=zw2(ig,2)
     638    endif
     639enddo
     640!
     641
     642!==============================================================================
     643!boucle de calcul de la vitesse verticale dans le thermique
     644!==============================================================================
     645do l=2,klev-1
     646!==============================================================================
     647
     648
     649! On decide si le thermique est encore actif ou non
     650! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     651    do ig=1,ngrid
     652       active(ig)=active(ig) &
     653&                 .and. zw2(ig,l)>1.e-10 &
     654&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     655    enddo
     656
     657
     658
     659!---------------------------------------------------------------------------
     660! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     661! sans tenir compte du detrainement et de l'entrainement dans cette
     662! couche
     663! C'est a dire qu'on suppose
     664! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
     665! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     666! avant) a l'alimentation pour avoir un calcul plus propre
     667!---------------------------------------------------------------------------
     668
     669   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
     670   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
     671
     672    do ig=1,ngrid
     673!       print*,'active',active(ig),ig,l
     674        if(active(ig)) then
     675        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
     676        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
     677        zta_est(ig,l)=ztva_est(ig,l)
     678        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
     679        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
     680     &      -zqla_est(ig,l))-zqla_est(ig,l))
     681
     682!------------------------------------------------
     683!AJAM:nouveau calcul de w² 
     684!------------------------------------------------
     685              zdz=zlev(ig,l+1)-zlev(ig,l)
     686              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     687
     688              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     689              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
     690              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
     691 
     692
     693             if (w_est(ig,l+1).lt.0.) then
     694                w_est(ig,l+1)=zw2(ig,l)
     695             endif
     696       endif
     697    enddo
     698
     699
     700!-------------------------------------------------
     701!calcul des taux d'entrainement et de detrainement
     702!-------------------------------------------------
     703
     704     do ig=1,ngrid
     705        if (active(ig)) then
     706
     707          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
     708          zw2m=w_est(ig,l+1)
     709          zdz=zlev(ig,l+1)-zlev(ig,l)
     710          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     711!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
     712          zbuoybis=zbuoy(ig,l)
     713          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
     714          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
     715
     716         
     717          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
     718    &     afact*zbuoybis/zw2m - fact_epsilon )
     719
     720
     721          detr_star(ig,l)=f_star(ig,l)*zdz                        &
     722    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
     723    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
     724         
     725! En dessous de lalim, on prend le max de alim_star et entr_star pour
     726! alim_star et 0 sinon
     727        if (l.lt.lalim(ig)) then
     728          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
     729          entr_star(ig,l)=0.
     730        endif
     731
     732! Calcul du flux montant normalise
     733      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
     734     &              -detr_star(ig,l)
     735
     736      endif
     737   enddo
     738
     739
     740!----------------------------------------------------------------------------
     741!calcul de la vitesse verticale en melangeant Tl et qt du thermique
     742!---------------------------------------------------------------------------
     743   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
     744   do ig=1,ngrid
     745       if (activetmp(ig)) then
     746           Zsat=.false.
     747           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
     748     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
     749     &            /(f_star(ig,l+1)+detr_star(ig,l))
     750           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
     751     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
     752     &            /(f_star(ig,l+1)+detr_star(ig,l))
     753
     754        endif
     755    enddo
     756
     757   ztemp(:)=zpspsk(:,l)*ztla(:,l)
     758   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
     759
     760   do ig=1,ngrid
     761      if (activetmp(ig)) then
     762        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
     763! on ecrit de maniere conservative (sat ou non)
     764!          T = Tl +Lv/Cp ql
     765           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
     766           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
     767           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
     768!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
     769           zha(ig,l) = ztva(ig,l)
     770           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
     771     &              -zqla(ig,l))-zqla(ig,l))
     772           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     773           zdz=zlev(ig,l+1)-zlev(ig,l)
     774           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
     775
     776            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     777            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
     778            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
     779      endif
     780   enddo
     781
     782        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
     783!
     784!---------------------------------------------------------------------------
     785!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
     786!---------------------------------------------------------------------------
     787
     788   do ig=1,ngrid
     789            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
     790!               stop'On tombe sur le cas particulier de thermcell_dry'
     791                print*,'On tombe sur le cas particulier de thermcell_plume'
     792                zw2(ig,l+1)=0.
     793                linter(ig)=l+1
     794            endif
     795
     796        if (zw2(ig,l+1).lt.0.) then
     797           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
     798     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
     799           zw2(ig,l+1)=0.
     800        endif
     801
     802           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
     803
     804        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
     805!   lmix est le niveau de la couche ou w (wa_moy) est maximum
     806!on rajoute le calcul de lmix_bis
     807            if (zqla(ig,l).lt.1.e-10) then
     808               lmix_bis(ig)=l+1
     809            endif
     810            lmix(ig)=l+1
     811            wmaxa(ig)=wa_moy(ig,l+1)
     812        endif
     813   enddo
     814
     815!=========================================================================
     816! FIN DE LA BOUCLE VERTICALE
     817      enddo
     818!=========================================================================
     819
     820!on recalcule alim_star_tot
     821       do ig=1,ngrid
     822          alim_star_tot(ig)=0.
     823       enddo
     824       do ig=1,ngrid
     825          do l=1,lalim(ig)-1
     826          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     827          enddo
     828       enddo
     829       
     830
     831        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
     832
     833#undef wrgrads_thermcell
     834#ifdef wrgrads_thermcell
     835         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
     836         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
     837         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
     838         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
     839         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
     840         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
     841         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
     842#endif
     843
     844
     845     return
     846     end
     847
Note: See TracChangeset for help on using the changeset viewer.