Changeset 1399 for LMDZ4/branches


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
Location:
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

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

    r1311 r1399  
    88     &      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,  &
    99     &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
    10      &       zmax0,f0,zw2,fraca)
     10     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl)
    1111
    1212      USE dimphy
    … …  
    4949      real zqla(klon,klev)
    5050      real zqta(klon,klev)
     51      real ztv(klon,klev)
     52      real zpspsk(klon,klev)
     53      real ztla(klon,klev)
     54      real zthl(klon,klev)
    5155      real wmax_sec(klon)
    5256      real zmax_sec(klon)
    … …  
    214218     &      ,r_aspect_thermals,l_mix_thermals &
    215219     &      ,tau_thermals,iflag_thermals_ed,Ale,Alp,lalim_conv,wght_th &
    216      &      ,zmax0,f0,zw2,fraca)
     220     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
     221     &      ,ztla,zthl)
    217222           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
    218223         else

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F

    r1299 r1399  
    77     s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
    88     s                   frac_impa, frac_nucl,
    9      s                   prfl, psfl, rhcl)
     9     s                   prfl, psfl, rhcl, zqta, fraca,
     10     s                   ztv, zpspsk, ztla, zthl, iflag_cldcon)
    1011
    1112c
    … …  
    4142      REAL snow(klon) ! neige (mm/s)
    4243      REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    43       REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     44      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     45      REAL ztv(klon,klev)
     46      REAL zqta(klon,klev),fraca(klon,klev)
     47      REAL sigma1(klon,klev),sigma2(klon,klev)
     48      REAL qltot(klon,klev),ctot(klon,klev)
     49      REAL zpspsk(klon,klev),ztla(klon,klev)
     50      REAL zthl(klon,klev)
     51
    4452cAA
    4553c Coeffients de fraction lessivee : pour OFF-LINE
    … …  
    6371
    6472      INTEGER ninter ! sous-intervals pour la precipitation
    65       INTEGER ncoreczq
     73      INTEGER ncoreczq 
     74      INTEGER iflag_cldcon
    6675      PARAMETER (ninter=5)
    6776      LOGICAL evap_prec ! evaporation de la pluie
    … …  
    7281      real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    7382      real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
    74       real erf
     83      real erf   
     84      REAL qcloud(klon)
    7585c
    7686      LOGICAL cpartiel ! condensation partielle
    … …  
    8292c
    8393      INTEGER i, k, n, kk
    84       REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
     94      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
    8595      REAL zrfl(klon), zrfln(klon), zqev, zqevt
    8696      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
    … …  
    130140      zdelq=0.0
    131141     
     142      print*,'CLOUDTH4 A. JAM'
    132143      IF (appel1er) THEN
    133144c
    … …  
    322333c de l'eau condensee:
    323334c
     335
    324336      IF (cpartiel) THEN
    325337
    … …  
    351363                zq(i)=1.e-15
    352364              endif
    353            enddo
    354            do i=1,klon
     365           enddo
     366
     367              if (iflag_cldcon.eq.5) then
     368
     369                 call cloudth(klon,klev,k,ztv,
     370     .           zq,zqta,fraca,
     371     .           qcloud,ctot,zpspsk,paprs,ztla,zthl,
     372     .           ratqs,zqs,t)
     373
     374                 do i=1,klon
     375                 rneb(i,k)=ctot(i,k)
     376                 zqn(i)=qcloud(i)
     377                 enddo
     378
     379              else
     380
     381            do i=1,klon
    355382            zpdf_sig(i)=ratqs(i,k)*zq(i)
    356383            zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
    … …  
    372399            endif
    373400           
    374            enddo
     401           enddo
     402
     403         endif ! iflag_cldcon
    375404
    376405        endif ! iflag_pdf

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F

    r1392 r1399  
    669669
    670670      REAL zw2(klon,klev+1)
    671       REAL fraca(klon,klev+1)
     671      REAL fraca(klon,klev+1)       
     672      REAL ztv(klon,klev)
     673      REAL zpspsk(klon,klev)
     674      REAL ztla(klon,klev)
     675      REAL zthl(klon,klev)
    672676
    673677c Variables locales pour la couche limite (al1):
    … …  
    24402444     s      ,ratqsdiff,zqsatth
    24412445con rajoute ale et alp, et les caracteristiques de la couche alim
    2442      s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
     2446     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
     2447     s      ,ztv,zpspsk,ztla,zthl)
    24432448
    24442449! ----------------------------------------------------------------------
    … …  
    27022707     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
    27032708     .           frac_impa, frac_nucl,
    2704      .           prfl, psfl, rhcl)
     2709     .           prfl, psfl, rhcl,
     2710     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
    27052711
    27062712      WHERE (rain_lsc < 0) rain_lsc = 0.

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

    r1377 r1399  
    1010     &                  ,r_aspect,l_mix,tau_thermals,iflag_thermals_ed &
    1111     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
    12      &                  ,zmax0, f0,zw2,fraca)
     12     &                  ,zmax0, f0,zw2,fraca,ztv &
     13     &                  ,zpspsk,ztla,zthl)
    1314
    1415      USE dimphy
    … …  
    389390!      print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
    390391      if (iflag_thermals_ed<=9) then
    391 !         print*,'THERM NOVUELLE/NOUVELLE/ANCIENNE'
     392!         print*,'THERM NOUVELLE/NOUVELLE Arnaud'
    392393         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
    393394     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    … …  
    396397     &    ,lev_out,lunout1,igout)
    397398
    398       elseif (iflag_thermals_ed<=19) then
     399      elseif (iflag_thermals_ed>9) then
    399400!        print*,'THERM RIO et al 2010, version d Arnaud'
    400401         CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&

  • 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.