Changeset 1968


Ignore:
Timestamp:
Feb 11, 2014, 4:01:12 PM (11 years ago)
Author:
fhourdin
Message:

Nouvelle version de thermcell_plume, permettant d'avoir des stratocumulus

avec 40 couches dans le cas fire.

New version of thermcell_plume, giving stratocumulus in the fire case

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1907 r1968  
    99
    1010!--------------------------------------------------------------------------
    11 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
     11! thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
     12! Last modified : Arnaud Jam 2014/02/11
     13!                 Better representation of stratocumulus
    1214!--------------------------------------------------------------------------
    1315
    14       IMPLICIT NONE
     16       IMPLICIT NONE
    1517
    1618#include "YOMCST.h"
    … …  
    3840      integer lev_out                           ! niveau pour les print
    3941      integer nbpb
    40       real zcon2(ngrid)
    4142   
    4243      real alim_star_tot(ngrid)
    … …  
    6263
    6364      REAL ztva_est(ngrid,klev)
     65      REAL ztv_est(ngrid,klev)
    6466      REAL zqla_est(ngrid,klev)
    6567      REAL zqsatth(ngrid,klev)
    6668      REAL zta_est(ngrid,klev)
    67       REAL zdw2
     69      REAL ztemp(ngrid),zqsat(ngrid)
     70      REAL zdw2,zdw2bis
    6871      REAL zw2modif
    69       REAL zeps
     72      REAL zw2fact,zw2factbis
     73      REAL zeps(ngrid,klev)
    7074
    7175      REAL linter(ngrid)
    … …  
    7478      REAL    wmaxa(ngrid)
    7579
    76       INTEGER ig,l,k
    77 
    78       real zdz,zfact,zbuoy,zalpha,zdrag
     80      INTEGER ig,l,k,lt,it
     81
     82      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
     83      real zbuoyjam(ngrid,klev)
     84      real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
    7985      real zcor,zdelta,zcvm5,qlbef
    80       real Tbef,qsatbef
    81       real dqsat_dT,DT,num,denom
     86      real betalpha,zbetalpha
     87      real eps, afact
    8288      REAL REPS,RLvCp,DDT0
    8389      PARAMETER (DDT0=.01)
    8490      logical Zsat
    8591      LOGICAL active(ngrid),activetmp(ngrid)
    86       REAL fact_gamma,fact_epsilon,fact_gamma2
     92      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
    8793      REAL c2(ngrid,klev)
    88       REAL a1,m
    89 
    90       REAL zw2fact,expa
    9194      Zsat=.false.
    9295! Initialisation
     96
    9397      RLvCp = RLVTT/RCPD
    94      
    95      
    96          fact_epsilon=0.002
    97          a1=2./3.
    98          fact_gamma=0.9
    99          zfact=fact_gamma/(1+fact_gamma)
    100          fact_gamma2=zfact
    101          expa=0.
    102      
    103 
    104 ! Initialisations des variables reeles
     98      fact_epsilon=0.002
     99      betalpha=0.9
     100      afact=2./3.           
     101
     102      zbetalpha=betalpha/(1.+betalpha)
     103
     104
     105! Initialisations des variables r?elles
    105106if (1==1) then
    106107      ztva(:,:)=ztv(:,:)
    107108      ztva_est(:,:)=ztva(:,:)
     109      ztv_est(:,:)=ztv(:,:)
    108110      ztla(:,:)=zthl(:,:)
    109111      zqta(:,:)=po(:,:)
     112      zqla(:,:)=0.
    110113      zha(:,:) = ztva(:,:)
    111114else
    112115      ztva(:,:)=0.
     116      ztv_est(:,:)=0.
    113117      ztva_est(:,:)=0.
    114118      ztla(:,:)=0.
    … …  
    128132      entr(:,:)=0.
    129133      zw2(:,:)=0.
     134      zbuoy(:,:)=0.
     135      zbuoyjam(:,:)=0.
     136      gamma(:,:)=0.
     137      zeps(:,:)=0.
    130138      w_est(:,:)=0.
    131139      f_star(:,:)=0.
    132140      wa_moy(:,:)=0.
    133141      linter(:)=1.
    134       linter(:)=1.
    135 
     142!     linter(:)=1.
    136143! Initialisation des variables entieres
    137144      lmix(:)=1
    … …  
    139146      wmaxa(:)=0.
    140147      lalim(:)=1
     148
    141149
    142150!-------------------------------------------------------------------------
    … …  
    156164               lalim(ig)=l+1
    157165               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     166               print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
    158167            endif
    159168         enddo
    … …  
    169178
    170179
     180
     181
    171182!------------------------------------------------------------------------------
    172183! Calcul dans la premiere couche
    173184! On decide dans cette version que le thermique n'est actif que si la premiere
    174185! couche est instable.
    175 ! Pourrait etre change si on veut que le thermiques puisse se déclencher
     186! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    176187! dans une couche l>1
    177188!------------------------------------------------------------------------------
    … …  
    210221
    211222
    212 ! Premier calcul de la vitesse verticale a partir de la temperature
    213 ! potentielle virtuelle
    214 !     if (1.eq.1) then
    215 !         w_est(ig,3)=zw2(ig,2)* &
    216 !    &      ((f_star(ig,2))**2) &
    217 !    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
    218 !    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
    219 !    &      *(zlev(ig,3)-zlev(ig,2))
    220 !      endif
    221 
    222 
    223223!---------------------------------------------------------------------------
    224224! calcul des proprietes thermodynamiques et de la vitesse de la couche l
    225225! sans tenir compte du detrainement et de l'entrainement dans cette
    226226! couche
     227! C'est a dire qu'on suppose
     228! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    227229! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    228230! avant) a l'alimentation pour avoir un calcul plus propre
    229231!---------------------------------------------------------------------------
    230232
    231      call thermcell_condens(ngrid,active, &
    232 &          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
    233 
    234     do ig=1,ngrid
    235         if(active(ig)) then
     233   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
     234   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
     235    do ig=1,ngrid
     236       print*,'active',active(ig),ig,l
     237        if(active(ig)) then
     238        zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
    236239        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    237240        zta_est(ig,l)=ztva_est(ig,l)
    … …  
    239242        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    240243     &      -zqla_est(ig,l))-zqla_est(ig,l))
    241 
    242          if (1.eq.0) then
    243 !calcul de w_est sans prendre en compte le drag
    244              w_est(ig,l+1)=zw2(ig,l)*  &
    245      &                   ((f_star(ig,l))**2)  &
    246      &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
    247      &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    248      &                   *(zlev(ig,l+1)-zlev(ig,l))
    249          else
    250 
    251            zdz=zlev(ig,l+1)-zlev(ig,l)
    252            zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
    253            zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    254            zdrag=fact_epsilon/(zalpha**expa)
    255            zw2fact=zbuoy/zdrag*a1
    256            w_est(ig,l+1)=(w_est(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
    257       &    +zw2fact
    258 
    259          endif
    260 
     244 
     245
     246!------------------------------------------------
     247!AJAM:nouveau calcul de w? 
     248!------------------------------------------------
     249              zdz=zlev(ig,l+1)-zlev(ig,l)
     250              zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
     251              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     252
     253              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     254              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
     255              zdw2=afact*zbuoy(ig,l)/fact_epsilon
     256              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
     257!             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
     258!             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
     259!    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     260!    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
     261!              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
     262             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     263    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     264    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
    261265             if (w_est(ig,l+1).lt.0.) then
    262                 w_est(ig,l+1)=zw2(ig,l)
     266!               w_est(ig,l+1)=zw2(ig,l)
     267                w_est(ig,l+1)=0.0001
    263268             endif
    264269       endif
    265270    enddo
    266271
     272
    267273!-------------------------------------------------
    268274!calcul des taux d'entrainement et de detrainement
    … …  
    271277     do ig=1,ngrid
    272278        if (active(ig)) then
     279
     280!          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
     281          zw2m=w_est(ig,l+1)
     282!          zw2m=zw2(ig,l)
    273283          zdz=zlev(ig,l+1)-zlev(ig,l)
    274           zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    275  
    276 ! estimation de la fraction couverte par les thermiques
    277           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
    278 
    279 !calcul de la soumission papier
    280 ! Calcul  du taux d'entrainement entr_star (epsilon)
    281            entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
    282      &     a1*zbuoy/w_est(ig,l+1) &
    283      &     - fact_epsilon/zalpha**expa  ) &
    284      &     +0. )
    285            
    286 !calcul du taux de detrainment (delta)
    287 !           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
    288 !     &      MAX(1.e-3, &
    289 !     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
    290 !     &      +0.01*(max(zqta(ig,l-1)-po(ig,l),0.)/(po(ig,l))/(w_est(ig,l+1)))**0.5    &   
    291 !     &     +0. ))
    292 
    293           m=0.5
    294 
    295           detr_star(ig,l)=1.*f_star(ig,l)*zdz *                    &
    296     &     MAX(5.e-4,-fact_gamma2*a1*(1./w_est(ig,l+1))*((1.-(1.-m)/(1.+70*zqta(ig,l-1)))*zbuoy        &
    297     &     -40*(1.-m)*(max(zqta(ig,l-1)-po(ig,l),0.))/(1.+70*zqta(ig,l-1)) )   )
    298 
    299 !           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
    300 !     &      MAX(0.0, &
    301 !     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
    302 !     &      +20*(max(zqta(ig,l-1)-po(ig,l),0.))**1*(zalpha/w_est(ig,l+1))**0.5    &   
    303 !     &     +0. ))
    304 
    305 
     284          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     285!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
     286          zbuoybis=zbuoy(ig,l)
     287          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
     288          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
     289
     290         
     291!          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
     292!    &     afact*zbuoybis/zw2m - fact_epsilon )
     293
     294!          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
     295!    &     afact*zbuoybis/zw2m - fact_epsilon )
     296
     297!Modif AJAM
     298         
     299        lmel=0.1*zlev(ig,l)
     300!        lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
     301        lt=l+1
     302         do it=1,klev-(l+1)
     303          zdz2=zlev(ig,lt)-zlev(ig,l)
     304          if (zdz2.gt.lmel) then
     305          zdz3=zlev(ig,lt)-zlev(ig,lt-1)
     306!          ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
     307!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
     308
     309         zbuoyjam(ig,l)=1.*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
     310    &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
     311    &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
     312
     313!         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
     314!    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
     315!    &          po(ig,lt-1))/po(ig,lt-1))
     316         
     317          else
     318          lt=lt+1
     319          endif
     320          enddo
     321
     322!          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     323
     324          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
     325    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
     326
     327          entrbis=entr_star(ig,l)
     328
     329
     330          detr_star(ig,l)=f_star(ig,l)*zdz                        &
     331    &     *MAX(1.e-4, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m          &
     332    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
     333
     334
     335!           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     336!
     337!           entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
     338!     &     afact*zbuoy(ig,l)/zw2m &
     339!     &     - 1.*fact_epsilon)
     340
     341         
    306342! En dessous de lalim, on prend le max de alim_star et entr_star pour
    307343! alim_star et 0 sinon
    … …  
    310346          entr_star(ig,l)=0.
    311347        endif
    312 
    313 !attention test
    314 !        if (detr_star(ig,l).gt.(f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l))) then       
    315 !            detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
     348!        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
     349!          alim_star(ig,l)=entrbis
    316350!        endif
     351
     352        print*,'alim0',l,lalim(ig),alim_star(ig,l),entrbis,f_star(ig,l)
    317353! Calcul du flux montant normalise
    318354      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    … …  
    321357      endif
    322358   enddo
     359
    323360
    324361!----------------------------------------------------------------------------
    … …  
    339376    enddo
    340377
    341    call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
    342 
    343 
     378   ztemp(:)=zpspsk(:,l)*ztla(:,l)
     379   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
    344380   do ig=1,ngrid
    345381      if (activetmp(ig)) then
    346382! on ecrit de maniere conservative (sat ou non)
    347383!          T = Tl +Lv/Cp ql
     384           zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
    348385           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    349386           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
    … …  
    352389           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
    353390     &              -zqla(ig,l))-zqla(ig,l))
    354 
    355 !on ecrit zqsat
    356            zqsatth(ig,l)=qsatbef 
    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) &
    363 !     &                 *(zlev(ig,l+1)-zlev(ig,l))
    364 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    365 ! La meme en plus modulaire :
    366            zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     391           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    367392           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 
    372 !if (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
    376 !else
    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 !endif
    382 
     393           zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
     394           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
     395
     396            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     397            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
     398            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
     399            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
     400!            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
     401            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     402    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     403    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
    383404      endif
    384405   enddo
    … …  
    441462        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    442463
    443 
    444       return
    445       end
     464#undef wrgrads_thermcell
     465#ifdef wrgrads_thermcell
     466         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
     467         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
     468         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
     469         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
     470         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
     471         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
     472         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
     473#endif
     474
     475
     476     return
     477     end
     478
    446479
    447480!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    625658! On decide dans cette version que le thermique n'est actif que si la premiere
    626659! couche est instable.
    627 ! Pourrait etre change si on veut que le thermiques puisse se déclencher
     660! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    628661! dans une couche l>1
    629662!------------------------------------------------------------------------------
    … …  
    686719
    687720!------------------------------------------------
    688 !AJAM:nouveau calcul de w² 
     721!AJAM:nouveau calcul de w? 
    689722!------------------------------------------------
    690723              zdz=zlev(ig,l+1)-zlev(ig,l)
    … …  
    855888     return
    856889     end
    857 
Note: See TracChangeset for help on using the changeset viewer.