Ignore:
Timestamp:
Jun 8, 2010, 10:29:11 AM (14 years ago)
Author:
idelkadi
Message:

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

  • calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

  • fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

  • physiq.F : appel modifie a calltherm et fisrtilp
  • thermcell_main.F90 : changement dans les flag_thermals_ed
  • thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

  • cloudth.F90 : PDF d'eau bi-gaussiennes
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1373 r1399  
    6464      REAL zqsatth(ngrid,klev)
    6565      REAL zta_est(ngrid,klev)
    66       REAL ztemp(ngrid),zqsat(ngrid)
    6766      REAL zdw2
    6867      REAL zw2modif
    … …  
    7675      INTEGER ig,l,k
    7776
    78       real zdz,zfact,zbuoy,zalpha
     77      real zdz,zfact,zbuoy,zalpha,zdrag
    7978      real zcor,zdelta,zcvm5,qlbef
     79      real Tbef,qsatbef
     80      real dqsat_dT,DT,num,denom
    8081      REAL REPS,RLvCp,DDT0
    8182      PARAMETER (DDT0=.01)
    … …  
    8485      REAL fact_gamma,fact_epsilon,fact_gamma2
    8586      REAL c2(ngrid,klev)
     87      REAL a1,m
    8688
    8789      REAL zw2fact,expa
    … …  
    9092      RLvCp = RLVTT/RCPD
    9193     
    92       if (iflag_thermals_ed==0) then
    93          fact_gamma=1.
    94          fact_epsilon=1.
    95       else if (iflag_thermals_ed==1)  then
    96 ! Valeurs au moment de la premiere soumission des papiers
    97          fact_gamma=1.
     94     
    9895         fact_epsilon=0.002
    99          fact_gamma2=0.6
    100 ! Suggestions des Fleurs, Septembre 2009
    101          fact_epsilon=0.015
    102 !test cr
    103 !         fact_epsilon=0.002
     96         a1=2./3.
    10497         fact_gamma=0.9
    105          fact_gamma2=0.7
    106 
    107       else if (iflag_thermals_ed==2)  then
    108          fact_gamma=1.
    109          fact_epsilon=2.
    110       else if (iflag_thermals_ed==3)  then
    111          fact_gamma=3./4.
    112          fact_epsilon=3.
    113       endif
    114 
    115 !      write(lunout,*)'THERM 31H '
    116 
    117       print*,'THERMCELL_PLUME OPTIMISE V0 '
     98         zfact=fact_gamma/(1+fact_gamma)
     99         fact_gamma2=zfact
     100         expa=0.
     101     
     102      print*,'THERM 31H '
     103
    118104! Initialisations des variables reeles
    119 if (1==0) then
     105if (1==1) then
    120106      ztva(:,:)=ztv(:,:)
    121107      ztva_est(:,:)=ztva(:,:)
    … …  
    168154               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
    169155     &                       *sqrt(zlev(ig,l+1))
    170                lalim(ig)=l+1
     156               lalim(:)=l+1
    171157               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    172158            endif
    … …  
    187173! On decide dans cette version que le thermique n'est actif que si la premiere
    188174! couche est instable.
    189 ! Pourrait etre change si on veut que le thermiques puisse se déclencher
     175! Pourrait etre change si on veut que le thermiques puisse se déclencher
    190176! dans une couche l>1
    191177!------------------------------------------------------------------------------
    … …  
    239225! sans tenir compte du detrainement et de l'entrainement dans cette
    240226! couche
    241 ! C'est a dire qu'on suppose
    242 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
    243227! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    244228! avant) a l'alimentation pour avoir un calcul plus propre
    245229!---------------------------------------------------------------------------
    246230
    247    ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
    248    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat)
    249 
     231     call thermcell_condens(ngrid,active, &
     232&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
    250233
    251234    do ig=1,ngrid
    252235        if(active(ig)) then
    253         zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsat(ig))
    254236        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
    255237        zta_est(ig,l)=ztva_est(ig,l)
    … …  
    258240     &      -zqla_est(ig,l))-zqla_est(ig,l))
    259241
     242         if (1.eq.0) then
     243!calcul de w_est sans prendre en compte le drag
    260244             w_est(ig,l+1)=zw2(ig,l)*  &
    261245     &                   ((f_star(ig,l))**2)  &
    … …  
    263247     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    264248     &                   *(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
    265261             if (w_est(ig,l+1).lt.0.) then
    266262                w_est(ig,l+1)=zw2(ig,l)
    … …  
    277273          zdz=zlev(ig,l+1)-zlev(ig,l)
    278274          zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    279           zfact=fact_gamma/(1.+fact_gamma)
    280275 
    281276! estimation de la fraction couverte par les thermiques
    282277          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
    283278
    284 !calcul de la soumission papier
    285           if (1.eq.1) then
    286              fact_epsilon=0.0007
    287 !             zfact=0.9/(1.+0.9)
    288              zfact=0.3
    289              fact_gamma=0.7
    290              fact_gamma2=0.6
    291              expa=0.25
     279!calcul de la soumission papier
    292280! Calcul  du taux d'entrainement entr_star (epsilon)
    293281           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
    294      &     zbuoy/w_est(ig,l+1) )&
    295 !- fact_epsilon/zalpha**0.25  ) &
    296      &     +0.000 )
    297 
    298 !           entr_star(ig,l)=f_star(ig,l)*zdz * (  1./3 * MAX(0.,  &     
    299 !     &     zbuoy/w_est(ig,l+1) - 1./zalpha**0.25  ) &
    300 !     &     +0.000 )
    301 ! Calcul du taux de detrainement detr_star (delta)
    302 !           if (zqla_est(ig,l).lt.1.e-10) then
     282     &     a1*zbuoy/w_est(ig,l+1) &
     283     &     - fact_epsilon/zalpha**expa  ) &
     284     &     +0. )
     285           
     286!calcul du taux de detrainment (delta)
    303287!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
    304 !     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
    305 !     &     +0.0006 )
    306 !           else
     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
    307299!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
    308 !     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
    309 !     &     +0.002 )
    310 !           endif
    311            detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
    312      &     fact_gamma2 * MAX(0., &
    313 !fact_epsilon/zalpha**0.25
    314      &      -zbuoy/w_est(ig,l+1) )       &
    315 !     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
    316 !test
    317      &     +  0.006*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &   
    318      &     +0.0000 )
    319           else
    320 
    321 ! Calcul  du taux d'entrainement entr_star (epsilon)
    322            entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
    323      &     zbuoy/w_est(ig,l+1) - fact_epsilon  ) &
    324      &     +0.0000 )
    325 
    326 ! Calcul du taux de detrainement detr_star (delta)
    327            detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
    328      &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
    329      &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
    330      &     +0.0000 )
    331 
    332           endif
    333 !endif choix du calcul de E* et D*
    334 
    335 !cr test
    336 !           entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l)     
    337 
    338 ! Prise en compte de la fraction
    339 !      detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5))
     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
    340305
    341306! En dessous de lalim, on prend le max de alim_star et entr_star pour
    … …  
    346311        endif
    347312
     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)
     316!        endif
    348317! Calcul du flux montant normalise
    349318      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    … …  
    370339    enddo
    371340
    372    ztemp(:)=zpspsk(:,l)*ztla(:,l)
    373    call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
     341   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
     342
    374343
    375344   do ig=1,ngrid
    … …  
    377346        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
    378347! on ecrit de maniere conservative (sat ou non)
    379            zqla(ig,l) = max(0.,zqta(ig,l)-zqsatth(ig,l))
    380348!          T = Tl +Lv/Cp ql
    381349           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
    … …  
    387355
    388356!on ecrit zqsat
     357           zqsatth(ig,l)=qsatbef 
    389358
    390359!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    407376           zw2(ig,l+1)=zw2modif+zdw2
    408377else
    409 ! Tentative de reecriture de l'equation de w2. A reprendre ...
    410 !           zdw2=2*zdz*zbuoy
    411 !           zw2modif=zw2(ig,l)*(1.-2.*zdz*(zeps+fact_epsilon))
    412 !!!!!      zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l))
    413 !formulation Arnaud
    414 !           zw2fact=zbuoy*zalpha**expa/fact_epsilon
    415 !           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa*(1+fact_gamma))*zdz) &
    416 !      &    +zw2fact
    417            if (zbuoy.gt.1.e-10) then
    418            zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma-zfact)
    419            zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
     378           zdrag=fact_epsilon/(zalpha**expa)
     379           zw2fact=zbuoy/zdrag*a1
     380           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
    420381      &    +zw2fact
    421            else
    422            zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma)
    423            zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
    424       &    +zw2fact
    425 
    426            endif
     382
    427383
    428384endif
    429 !           zw2(ig,l+1)=zw2modif+zdw2
     385
    430386      endif
    431387   enddo
    … …  
    440396            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    441397!               stop'On tombe sur le cas particulier de thermcell_dry'
    442                 write(lunout,*)                                         &
    443      &          'On tombe sur le cas particulier de thermcell_plume'
     398                print*,'On tombe sur le cas particulier de thermcell_plume'
    444399                zw2(ig,l+1)=0.
    445400                linter(ig)=l+1
    … …  
    487442      return
    488443      end
    489 
    490444
    491445!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracChangeset for help on using the changeset viewer.