Changeset 2140 for LMDZ5/trunk/libf


Ignore:
Timestamp:
Nov 7, 2014, 5:16:04 PM (10 years ago)
Author:
fhourdin
Message:

Nouvelle version des thermiques permettant de retrouver a 79 niveaux

des comportements proches de ceux obtenus dans la thèse d'Arnaud
Jam sur les cas de strato-cumulus.
thermcell_plume et cloudth modifiés.

New versions of thermal plume model for stratocumulus cases.

Arnaud Jam

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

Legend:

Unmodified
Added
Removed

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

    r1907 r2140  
    7474     
    7575!-----------------------------------------------------------------------------------------------------------------
    76 ! Initialisation des variables réelles
     76! Initialisation des variables r?elles
    7777!-----------------------------------------------------------------------------------------------------------------
    7878      sigma1(:,:)=0.
    … …  
    9797
    9898!------------------------------------------------------------------------------------------------------------------
    99 ! Calcul de la fraction du thermique et des écart-types des distributions
     99! Calcul de la fraction du thermique et des ?cart-types des distributions
    100100!------------------------------------------------------------------------------------------------------------------                 
    101101      do ind1=1,ngrid
    … …  
    140140
    141141!-----------------------------------------------------------------------------------------------------------------
    142 ! Calcul des écart-types pour s
    143 !-----------------------------------------------------------------------------------------------------------------
    144 
    145       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    146       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) 
    147 
    148  
    149 !-----------------------------------------------------------------------------------------------------------------
    150 ! Calcul de l'eau condensée et de la couverture nuageuse
     142! Calcul des ?cart-types pour s
     143!-----------------------------------------------------------------------------------------------------------------
     144
     145!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     146!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
     147!       if (paprs(ind1,ind2).gt.90000) then
     148!       ratqs(ind1,ind2)=0.002
     149!       else
     150!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
     151!       endif
     152       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
     153       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
     154!       sigma1s=ratqs(ind1,ind2)*po(ind1)
     155!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
     156 
     157!-----------------------------------------------------------------------------------------------------------------
     158! Calcul de l'eau condens?e et de la couverture nuageuse
    151159!-----------------------------------------------------------------------------------------------------------------
    152160      sqrt2pi=sqrt(2.*pi)
    … …  
    241249
    242250
    243 

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

    r2106 r2140  
    77     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    88     &           ,lev_out,lunout1,igout)
     9!    &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    910!--------------------------------------------------------------------------
    1011!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    … …  
    7677      REAL    wmaxa(ngrid)
    7778
    78       INTEGER ig,l,k,lt,it
     79      INTEGER ig,l,k,lt,it,lm
    7980
    8081      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
    … …  
    9798      REAL, SAVE :: betalpha, betalpha_omp=0.9
    9899      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
     100      REAL, SAVE :: fact_shell, fact_shell_omp=1.
     101      REAL,SAVE :: detr_min,detr_min_omp=1.e-5
     102      REAL,SAVE :: entr_min,entr_min_omp=1.e-5
    102103      REAL,SAVE :: detr_q_coef,detr_q_coef_omp=0.012
    103104      REAL,SAVE :: detr_q_power,detr_q_power_omp=0.5
    … …  
    286287 
    287288
    288 !------------------------------------------------
    289 !AJAM:nouveau calcul de w? 
    290 !------------------------------------------------
    291               zdz=zlev(ig,l+1)-zlev(ig,l)
    292               zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
    293               zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    294 
    295               zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    296               zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    297               zdw2=afact*zbuoy(ig,l)/fact_epsilon
    298               zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
    299 !             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
    300 !             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
    301 !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    302 !    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
    303 !              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
    304 
    305 !--------------------------------------------------
    306 !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
    307 !--------------------------------------------------
    308          if (iflag_thermals_ed==8) then
    309 ! Ancienne version
    310              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    311     &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    312     &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
    313  
    314 ! Nouvelle version Arnaud
    315          else
    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 
    326          endif
    327 
    328 
    329          if (iflag_thermals_ed<6) then
    330              zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    331 !              fact_epsilon=0.0005/(zalpha+0.025)**0.5
    332 !              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
    333               fact_epsilon=0.0002/(zalpha+0.1)
    334               zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
    335               zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
    336               zdw2=afact*zbuoy(ig,l)/fact_epsilon
    337               zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
    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))
    343 
    344          endif
    345 !--------------------------------------------------
    346 !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
    347 !on fait max(0.0001,.....)
    348 !--------------------------------------------------         
    349 
    350 !             if (w_est(ig,l+1).lt.0.) then
    351 !               w_est(ig,l+1)=zw2(ig,l)
    352 !                w_est(ig,l+1)=0.0001
    353 !             endif
    354 
    355        endif
    356     enddo
    357 
    358 
    359 !-------------------------------------------------
    360 !calcul des taux d'entrainement et de detrainement
    361 !-------------------------------------------------
    362 
    363      do ig=1,ngrid
    364         if (active(ig)) then
    365 
    366 !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
    367           zw2m=w_est(ig,l+1)
    368 !          zw2m=zw2(ig,l)
    369           zdz=zlev(ig,l+1)-zlev(ig,l)
    370           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    371 !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
    372           zbuoybis=zbuoy(ig,l)
    373           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
    374           zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
    375 
    376          
    377 !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
    378 !    &     afact*zbuoybis/zw2m - fact_epsilon )
    379 
    380 !          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
    381 !    &     afact*zbuoybis/zw2m - fact_epsilon )
    382 
    383289!Modif AJAM
    384          
     290
     291        zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     292        zdz=zlev(ig,l+1)-zlev(ig,l)         
    385293        lmel=fact_thermals_ed_dz*zlev(ig,l)
    386294!        lmel=0.09*zlev(ig,l)
    … …  
    451359                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
    452360    &                    +(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)
     361         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), &
     362   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    455363                elseif (zlmelup.ge.zinv) then
    456364                 ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
    … …  
    462370    &            ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
    463371
    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)
     372         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), &
     373   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    466374                else
    467375                   ztv_est(ig,l)=atv1*zlmel+btv1
    468376                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
    469377    &                           +(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)
     378                 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), &
     379   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    472380                endif
    473381
    … …  
    483391
    484392                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)
     393
     394                 print*,'on est pass?? par l??4',l,lt,ztv1,ztv2,ztv(ig,lt),ztv(ig,l),ztva_est(ig,l), &
     395   &                                   zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    487396!          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
    488397!    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
    … …  
    508417    &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
    509418    &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
    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)
     419         print*,'on est pass?? par l??',l,lt,ztv(ig,lt),ztva_est(ig,l),ztv(ig,l), &
     420   &                                   zbuoy(ig,l),zbuoyjam(ig,l)
    512421        endif !   if (iflag_thermals_ed.lt.8) then
     422
     423!------------------------------------------------
     424!AJAM:nouveau calcul de w? 
     425!------------------------------------------------
     426              zdz=zlev(ig,l+1)-zlev(ig,l)
     427              zdzbis=zlev(ig,l)-zlev(ig,l-1)
     428              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     429
     430              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     431              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
     432              zdw2=afact*zbuoy(ig,l)/fact_epsilon
     433              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
     434              lm=Max(1,l-2)
     435!              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
     436!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
     437!              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
     438!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
     439!             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
     440!             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
     441!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     442!    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
     443!              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
     444
     445!--------------------------------------------------
     446!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
     447!--------------------------------------------------
     448         if (iflag_thermals_ed==8) then
     449! Ancienne version
     450!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     451!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     452!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
     453
     454            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
     455
     456! Nouvelle version Arnaud
     457         else
     458!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     459!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     460!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
     461
     462            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
     463
     464!             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
     465!    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
     466!    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
     467
     468
     469
     470!            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
     471
     472!             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
     473!    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
     474
     475!             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
     476
     477         endif
     478
     479
     480         if (iflag_thermals_ed<6) then
     481             zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
     482!              fact_epsilon=0.0005/(zalpha+0.025)**0.5
     483!              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
     484              fact_epsilon=0.0002/(zalpha+0.1)
     485              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     486              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
     487              zdw2=afact*zbuoy(ig,l)/fact_epsilon
     488              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
     489!              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     490
     491!              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     492!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     493!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     494
     495            w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
     496
     497
     498         endif
     499!--------------------------------------------------
     500!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
     501!on fait max(0.0001,.....)
     502!--------------------------------------------------         
     503
     504!             if (w_est(ig,l+1).lt.0.) then
     505!               w_est(ig,l+1)=zw2(ig,l)
     506!                w_est(ig,l+1)=0.0001
     507!             endif
     508
     509       endif
     510    enddo
     511
     512
     513!-------------------------------------------------
     514!calcul des taux d'entrainement et de detrainement
     515!-------------------------------------------------
     516
     517     do ig=1,ngrid
     518        if (active(ig)) then
     519
     520!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
     521          zw2m=w_est(ig,l+1)
     522!          zw2m=zw2(ig,l)
     523          zdz=zlev(ig,l+1)-zlev(ig,l)
     524          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     525!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
     526          zbuoybis=zbuoy(ig,l)
     527          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
     528          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
     529
     530         
     531!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
     532!    &     afact*zbuoybis/zw2m - fact_epsilon )
     533
     534!          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
     535!    &     afact*zbuoybis/zw2m - fact_epsilon )
     536
     537
    513538
    514539!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    … …  
    526551          endif
    527552         
    528 !          if (zw2m.lt.1.) then
    529 !          zbuoyjam(ig,l)=zbuoy(ig,l)
    530 !          endif
    531 
    532553
    533554
    … …  
    537558    &     + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
    538559
     560!         detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
     561!    &                    ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
     562
    539563          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    540564
    … …  
    542566    &       mix0 * 0.1 / (zalpha+0.001)               &
    543567    &     + zbetalpha*MAX(entr_min,                   &
    544     &     afact*zbuoy(ig,l)/zw2m - fact_epsilon))
     568    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon ))
     569
     570
     571!          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
     572!    &       mix0 * 0.1 / (zalpha+0.001)               &
     573!    &     + MAX(entr_min,                   &
     574!    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  &
     575!    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
     576
     577
     578!         entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
     579!    &                    ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
    545580
    546581!          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
    … …  
    601636           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    602637           zdz=zlev(ig,l+1)-zlev(ig,l)
    603            zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
     638           zdzbis=zlev(ig,l)-zlev(ig,l-1)
    604639           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
    605640!!!!!!!          fact_epsilon=0.002
    … …  
    608643            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
    609644            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
     645!              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
     646!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
     647!              lm=Max(1,l-2)
     648!              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
     649!    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
    610650!            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
    611651!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
    612652!     &                   (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))
     653!            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
    614654!            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
    615655!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     656!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
     657            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
     658!            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
     659!    &                     (zw2(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdz+zdzbis))* &
    616660!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
     661
    617662
    618663           if (iflag_thermals_ed.lt.6) then
    … …  
    629674!    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
    630675!    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    631             zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     676!            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
     677            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
    632678
    633679           endif
    … …  
    759805&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
    760806&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    761 &           ,lev_out,lunout1,igout)
     807&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    762808
    763809!--------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.