Changeset 2106 for LMDZ5/trunk/libf


Ignore:
Timestamp:
Aug 11, 2014, 10:56:54 AM (10 years ago)
Author:
fhourdin
Message:

Version moidifiÃe des thermiques pour une meilleure représnation
des stratocumulus + contrôle par lecture de .def d'un ertain nombre
de paramètres internes.
Modified version of thermals for stratocumulus + reading new free parameters.
Modification by Arnaud Jam and Frederic Hourdin

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2046 r2106  
    1010!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    1111!--------------------------------------------------------------------------
     12USE IOIPSL, ONLY : getin
    1213
    1314       IMPLICIT NONE
    … …  
    7879
    7980      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
    80       real zbuoyjam(ngrid,klev)
     81      real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev)
    8182      real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
    82       real fact_shell
    8383      real ztv1,ztv2,factinv,zinv,zlmel
     84      real zlmelup,zlmeldwn,zlt,zltdwn,zltup
     85      real atv1,atv2,btv1,btv2
    8486      real ztv_est1,ztv_est2
    8587      real zcor,zdelta,zcvm5,qlbef
    86       real betalpha,zbetalpha
    87       real eps, afact
     88      real zbetalpha
     89      real eps
    8890      REAL REPS,RLvCp,DDT0
    8991      PARAMETER (DDT0=.01)
    9092      logical Zsat
    9193      LOGICAL active(ngrid),activetmp(ngrid)
    92       REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
     94      REAL fact_gamma,fact_gamma2,fact_epsilon2
     95
     96      REAL, SAVE :: fact_epsilon, fact_epsilon_omp=0.002
     97      REAL, SAVE :: betalpha, betalpha_omp=0.9
     98      REAL, SAVE :: afact, afact_omp=2./3.
     99      REAL, SAVE :: fact_shell, fact_shell_omp=0.6
     100      REAL,SAVE :: detr_min,detr_min_omp=1.e-4
     101      REAL,SAVE :: entr_min,entr_min_omp=1.e-4
     102      REAL,SAVE :: detr_q_coef,detr_q_coef_omp=0.012
     103      REAL,SAVE :: detr_q_power,detr_q_power_omp=0.5
     104      REAL,SAVE :: mix0,mix0_omp=0.
     105
     106      LOGICAL, SAVE :: first=.true.
     107
    93108      REAL c2(ngrid,klev)
     109
     110      if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
    94111      Zsat=.false.
    95112! Initialisation
    96113
    97114      RLvCp = RLVTT/RCPD
    98       fact_epsilon=0.002
    99       betalpha=0.9
    100       afact=2./3.   
    101       fact_shell=0.85       
     115      IF (first) THEN
     116     !$OMP MASTER
     117! FH : if ok_sync=.true. , the time axis is written at each time step
     118! in the output files. Only at the end in the opposite case
     119     CALL getin('thermals_fact_epsilon',fact_epsilon_omp)
     120     CALL getin('thermals_betalpha',betalpha_omp)
     121     CALL getin('thermals_afact',afact_omp)
     122     CALL getin('thermals_fact_shell',fact_shell_omp)
     123     CALL getin('thermals_detr_min',detr_min_omp)
     124     CALL getin('thermals_entr_min',entr_min_omp)
     125     CALL getin('thermals_detr_q_coef',detr_q_coef_omp)
     126     CALL getin('thermals_detr_q_power',detr_q_power_omp)
     127     CALL getin('thermals_mix0',mix0_omp)
     128!    CALL getin('thermals_X',X_omp)
     129!    X=X_omp
     130     !$OMP END MASTER
     131     !$OMP BARRIER
     132     fact_epsilon=fact_epsilon_omp
     133     betalpha=betalpha_omp
     134     afact=afact_omp
     135     fact_shell=fact_shell_omp
     136     detr_min=detr_min_omp
     137     entr_min=entr_min_omp
     138     detr_q_coef=detr_q_coef_omp
     139     detr_q_power=detr_q_power_omp
     140     mix0=mix0_omp
     141      first=.false.
     142      ENDIF
    102143
    103144      zbetalpha=betalpha/(1.+betalpha)
    … …  
    273314! Nouvelle version Arnaud
    274315         else
    275               w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    276     &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    277     &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     316!              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     317!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     318!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
     319
     320             w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     321
     322!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
     323!    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
     324!             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
     325
    278326         endif
    279327
    … …  
    288336              zdw2=afact*zbuoy(ig,l)/fact_epsilon
    289337              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
    290               w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    291     &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    292     &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     338              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     339
     340!              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     341!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     342!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    293343
    294344         endif
    … …  
    296346!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
    297347!on fait max(0.0001,.....)
    298 !--------------------------------------------------     
     348!--------------------------------------------------       
    299349
    300350!             if (w_est(ig,l+1).lt.0.) then
    … …  
    333383!Modif AJAM
    334384         
    335         lmel=fact_thermals_ed_dz*zlev(ig,l)
     385        lmel=fact_thermals_ed_dz*zlev(ig,l)
     386!        lmel=0.09*zlev(ig,l)
    336387        zlmel=zlev(ig,l)+lmel
    337 !        lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
     388        zlmelup=zlmel+(zdz/2)
     389        zlmeldwn=zlmel-(zdz/2)
     390
    338391        lt=l+1
    339         zdz2=zlev(ig,lt)-zlev(ig,l)
     392        zlt=zlev(ig,lt)
     393        zdz3=zlev(ig,lt+1)-zlt
     394        zltdwn=zlt-zdz3/2
     395        zltup=zlt+zdz3/2
     396         
     397!=========================================================================
     398! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
     399!=========================================================================
     400
    340401!--------------------------------------------------
    341 !AJ052014: J'ai remplac? la boucle do par un do while
     402        if (iflag_thermals_ed.lt.8) then
     403!--------------------------------------------------
     404!AJ052014: J'ai remplac?? la boucle do par un do while
    342405! afin de faire moins de calcul dans la boucle
    343406!--------------------------------------------------
    344          do while (zdz2.lt.lmel)
    345          lt=lt+1
    346          zdz2=zlev(ig,lt)-zlev(ig,l)
    347          end do
    348 
    349          zdz3=zlev(ig,lt)-zlev(ig,lt-1)
    350 
     407            do while (zlmelup.gt.zltup)
     408               lt=lt+1
     409               zlt=zlev(ig,lt)
     410               zdz3=zlev(ig,lt+1)-zlt
     411               zltdwn=zlt-zdz3/2
     412               zltup=zlt+zdz3/2       
     413            enddo
    351414!--------------------------------------------------
    352415!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
    353 ! on cherche o? se trouve l'altitude d'inversion
     416! on cherche o?? se trouve l'altitude d'inversion
    354417! en calculant ztv1 (interpolation de la valeur de
    355418! theta au niveau lt en utilisant les niveaux lt-1 et
    356419! lt-2) et ztv2 (interpolation avec les niveaux lt+1
    357 ! et lt+2). Si theta r?ellement calcul?e au niveau lt
     420! et lt+2). Si theta r??ellement calcul??e au niveau lt
    358421! comprise entre ztv1 et ztv2, alors il y a inversion
    359422! et on calcule son altitude zinv en supposant que ztv(lt)
    360 ! est une combinaison lin?aire de ztv1 et ztv2.
    361 ! Ensuite, on calcule la flottabilit? en comparant
    362 ! la temp?rature de la couche l ? celle de l'air situ?
    363 ! l+lmel plus haut, ce qui n?cessite de savoir quel fraction
     423! est une combinaison lineaire de ztv1 et ztv2.
     424! Ensuite, on calcule la flottabilite en comparant
     425! la temperature de la couche l a celle de l'air situe
     426! l+lmel plus haut, ce qui necessite de savoir quel fraction
    364427! de cet air est au-dessus ou en-dessous de l'inversion   
    365428!--------------------------------------------------
    366 
    367 
    368          if (iflag_thermals_ed.lt.8) then
    369 
    370           ztv1=(ztv(ig,lt-1)-ztv(ig,lt-2))*zlev(ig,lt)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
    371     &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
     429            atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
     430            btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
    372431    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
    373 
    374           ztv2=(ztv(ig,lt+2)-ztv(ig,lt+1))*zlev(ig,lt)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
    375     &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
     432            atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
     433            btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
    376434    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
    377435
    378           if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then 
    379 
    380           factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
    381           zinv=zlev(ig,lt-1)+factinv*(zlev(ig,lt)-zlev(ig,lt-1))
    382 
    383           if (zlmel+0.5*zdz.ge.zinv) then
    384            if (zlmel-0.5*zdz.ge.zinv) then
    385 
    386           ztv_est(ig,l)=(ztv(ig,lt+2)-ztv(ig,lt+1))*(zlmel-0.*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
    387     &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
    388     &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
    389 
    390           zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
    391 
    392           else
    393 
    394           ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
    395     &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
    396     &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
    397           ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
    398     &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
    399     &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
    400           zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- &
    401     &          ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- &
    402     &          ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
    403 
    404            endif
    405 
    406           else
    407 
    408           ztv_est(ig,l)=(ztv(ig,lt-1)-ztv(ig,lt-2))*(zlmel-0.*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
    409     &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
    410     &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
    411 
    412           zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
    413 !         ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
    414 !    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
    415 !    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
    416 !          ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
    417 !    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
    418 !    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
    419 !          zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- &
    420 !    &          ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- &
    421 !    &          ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
    422 
    423 
    424 
    425 
    426 !          print*,'on est pass? par l?',l,lt,zbuoyjam(ig,l),zbuoy(ig,l)
    427           endif
    428 
    429 
    430           else
     436             ztv1=atv1*zlt+btv1
     437             ztv2=atv2*zlt+btv2
     438
     439             if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then 
     440
     441!--------------------------------------------------
     442!AJ052014: D??calage de zinv qui est entre le haut
     443!          et le bas de la couche lt
     444!--------------------------------------------------
     445                factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
     446                zinv=zltdwn+zdz3*factinv
     447
    431448         
    432 !          ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
    433 !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
    434 
    435          zbuoyjam(ig,l)=fact_shell*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
    436     &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
    437     &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
    438 
     449                if (zlmeldwn.ge.zinv) then
     450                   ztv_est(ig,l)=atv2*zlmel+btv2
     451                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
     452    &                    +(1.-fact_shell)*zbuoy(ig,l)
     453!         print*,'on est pass?? par l??1',l,lt,ztv1,ztv2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
     454!   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
     455                elseif (zlmelup.ge.zinv) then
     456                 ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
     457                   ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1
     458                   ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1
     459
     460                   zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
     461    &            ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
     462    &            ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
     463
     464!         print*,'on est pass?? par l??2',l,lt,ztv_est1,ztv_est2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
     465!   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
     466                else
     467                   ztv_est(ig,l)=atv1*zlmel+btv1
     468                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
     469    &                           +(1.-fact_shell)*zbuoy(ig,l)
     470!                 print*,'on est pass?? par l??3',l,lt,ztv1,ztv2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
     471!   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
     472                endif
     473
     474             else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
     475
     476                if (zlmeldwn.gt.zltdwn) then
     477                   zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- &
     478    &                ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
     479                else
     480                   zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
     481    &                ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
     482    &                ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
     483
     484                endif
     485!                 print*,'on est pass?? par l??4',l,lt,ztv1,ztv2,ztv(ig,lt),ztv(ig,l),ztva_est(ig,l), &
     486!   &                                   zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
     487!          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
     488!    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
     489!    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
    439490!         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
    440491!    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
    441492!     &          po(ig,lt-1))/po(ig,lt-1))
    442 
    443           endif
    444 
    445           else
    446 
    447          zbuoyjam(ig,l)=1.*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
     493          endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
     494
     495        else  !   if (iflag_thermals_ed.lt.8) then
     496           lt=l+1
     497           zdz2=zlev(ig,lt)-zlev(ig,l)
     498
     499           do while (lmel.gt.zdz2)
     500             lt=lt+1
     501             zlt=zlev(ig,lt)
     502             zdz2=zlt-zlev(ig,l)
     503           enddo
     504           zdz3=zlev(ig,lt+1)-zlt
     505           zltdwn=zlev(ig,lt)-zdz3/2
     506
     507           zbuoyjam(ig,l)=1.*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
    448508    &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
    449509    &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
    450 
    451           endif
    452          
     510!         print*,'on est pass?? par l??',l,lt,ztv(ig,lt),ztva_est(ig,l),ztv(ig,l), &
     511!   &                                   zbuoy(ig,l),zbuoyjam(ig,l)
     512        endif !   if (iflag_thermals_ed.lt.8) then
    453513
    454514!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     515
     516!=========================================================================
     517! 4. Calcul de l'entrainement et du detrainement
     518!=========================================================================
    455519
    456520!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
    457521!    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
    458 
    459522!          entrbis=entr_star(ig,l)
    460523
    461           if (iflag_thermals_ed.lt.6) then
     524          if (iflag_thermals_ed.lt.6) then
    462525          fact_epsilon=0.0002/(zalpha+0.1)
    463           endif
     526          endif
     527         
     528!          if (zw2m.lt.1.) then
     529!          zbuoyjam(ig,l)=zbuoy(ig,l)
     530!          endif
     531
     532
    464533
    465534          detr_star(ig,l)=f_star(ig,l)*zdz             &
    466     &     *MAX(1.e-4, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m          &
    467     &     + 0.012*(zdqt(ig,l)/zw2m)**0.5)
     535    &     *( mix0 * 0.1 / (zalpha+0.001)               &
     536    &     + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m   &
     537    &     + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
    468538
    469539          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    470540
    471           entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
    472     &     afact*zbuoy(ig,l)/zw2m - fact_epsilon)
     541          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
     542    &       mix0 * 0.1 / (zalpha+0.001)               &
     543    &     + zbetalpha*MAX(entr_min,                   &
     544    &     afact*zbuoy(ig,l)/zw2m - fact_epsilon))
    473545
    474546!          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
    … …  
    487559!        endif
    488560
    489 !       print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
     561!        print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
    490562! Calcul du flux montant normalise
    491563      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    … …  
    496568
    497569
    498 !----------------------------------------------------------------------------
    499 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
    500 !---------------------------------------------------------------------------
     570!============================================================================
     571! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
     572!===========================================================================
     573
    501574   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
    502575   do ig=1,ngrid
    … …  
    530603           zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
    531604           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
    532            fact_epsilon=0.002
     605!!!!!!!          fact_epsilon=0.002
    533606            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    534607            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    535608            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
    536609            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
    537 !            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
    538             zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    539     &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    540     &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    541 
    542            if (iflag_thermals_ed.lt.6) then
     610!            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
     611!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
     612!     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
     613            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     614!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     615!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     616!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
     617
     618           if (iflag_thermals_ed.lt.6) then
    543619           zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
    544620!           fact_epsilon=0.0005/(zalpha+0.025)**0.5
    … …  
    550626            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
    551627
    552             zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    553     &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    554     &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    555 
    556            endif
     628!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     629!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     630!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     631            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     632
     633           endif
    557634
    558635
    … …  
    562639        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
    563640!
    564 !---------------------------------------------------------------------------
    565 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    566 !---------------------------------------------------------------------------
     641!===========================================================================
     642! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
     643!===========================================================================
    567644
    568645   nbpb=0
    … …  
    10331110     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    10341111           zw2(ig,l+1)=0.
     1112        elseif (f_star(ig,l+1).lt.0.) then
     1113           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
     1114     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
     1115           print*,"linter plume", linter(ig)
     1116           zw2(ig,l+1)=0.
    10351117        endif
    10361118
Note: See TracChangeset for help on using the changeset viewer.