Ignore:
Timestamp:
Feb 18, 2010, 2:14:02 PM (14 years ago)
Author:
Laurent Fairhead
Message:

Modifications to thermals for TKE transport

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1299 r1311  
    114114      endif
    115115
    116       write(lunout,*)'THERM 31H '
     116!      write(lunout,*)'THERM 31H '
    117117
    118118! Initialisations des variables reeles
    … …  
    483483      return
    484484      end
     485
     486
     487!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     488!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     489!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     490!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     491!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     492!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     493 SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
     494&           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
     495&           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     496&           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     497&           ,lev_out,lunout1,igout)
     498
     499!--------------------------------------------------------------------------
     500!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
     501! Version conforme a l'article de Rio et al. 2010.
     502! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
     503!--------------------------------------------------------------------------
     504
     505      IMPLICIT NONE
     506
     507#include "YOMCST.h"
     508#include "YOETHF.h"
     509#include "FCTTRE.h"
     510#include "iniprint.h"
     511#include "thermcell.h"
     512
     513      INTEGER itap
     514      INTEGER lunout1,igout
     515      INTEGER ngrid,klev
     516      REAL ptimestep
     517      REAL ztv(ngrid,klev)
     518      REAL zthl(ngrid,klev)
     519      REAL po(ngrid,klev)
     520      REAL zl(ngrid,klev)
     521      REAL rhobarz(ngrid,klev)
     522      REAL zlev(ngrid,klev+1)
     523      REAL pplev(ngrid,klev+1)
     524      REAL pphi(ngrid,klev)
     525      REAL zpspsk(ngrid,klev)
     526      REAL alim_star(ngrid,klev)
     527      REAL f0(ngrid)
     528      INTEGER lalim(ngrid)
     529      integer lev_out                           ! niveau pour les print
     530   
     531      real alim_star_tot(ngrid)
     532
     533      REAL ztva(ngrid,klev)
     534      REAL ztla(ngrid,klev)
     535      REAL zqla(ngrid,klev)
     536      REAL zqta(ngrid,klev)
     537      REAL zha(ngrid,klev)
     538
     539      REAL detr_star(ngrid,klev)
     540      REAL coefc
     541      REAL entr_star(ngrid,klev)
     542      REAL detr(ngrid,klev)
     543      REAL entr(ngrid,klev)
     544
     545      REAL csc(ngrid,klev)
     546
     547      REAL zw2(ngrid,klev+1)
     548      REAL w_est(ngrid,klev+1)
     549      REAL f_star(ngrid,klev+1)
     550      REAL wa_moy(ngrid,klev+1)
     551
     552      REAL ztva_est(ngrid,klev)
     553      REAL zqla_est(ngrid,klev)
     554      REAL zqsatth(ngrid,klev)
     555      REAL zta_est(ngrid,klev)
     556      REAL zdw2
     557      REAL zw2modif
     558      REAL zw2fact
     559      REAL zeps(ngrid,klev)
     560
     561      REAL linter(ngrid)
     562      INTEGER lmix(ngrid)
     563      INTEGER lmix_bis(ngrid)
     564      REAL    wmaxa(ngrid)
     565
     566      INTEGER ig,l,k
     567
     568      real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
     569      real zbuoybis
     570      real zcor,zdelta,zcvm5,qlbef,zdz2
     571      real betalpha,zbetalpha
     572      real Tbef,qsatbef,b1,eps, afact
     573      real dqsat_dT,DT,num,denom,m
     574      REAL REPS,RLvCp,DDT0
     575      PARAMETER (DDT0=.01)
     576      logical Zsat
     577      LOGICAL active(ngrid),activetmp(ngrid)
     578      REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
     579      REAL c2(ngrid,klev)
     580      Zsat=.false.
     581! Initialisation
     582
     583      RLvCp = RLVTT/RCPD
     584      fact_epsilon=0.002
     585      betalpha=0.9
     586      afact=2./3.           
     587
     588      zbetalpha=betalpha/(1.+betalpha)
     589
     590!      print*,'THERM 31B'
     591
     592! Initialisations des variables reeles
     593if (1==0) then
     594      ztva(:,:)=ztv(:,:)
     595      ztva_est(:,:)=ztva(:,:)
     596      ztla(:,:)=zthl(:,:)
     597      zqta(:,:)=po(:,:)
     598      zha(:,:) = ztva(:,:)
     599else
     600      ztva(:,:)=0.
     601      ztva_est(:,:)=0.
     602      ztla(:,:)=0.
     603      zqta(:,:)=0.
     604      zha(:,:) =0.
     605endif
     606
     607      zqla_est(:,:)=0.
     608      zqsatth(:,:)=0.
     609      zqla(:,:)=0.
     610      detr_star(:,:)=0.
     611      entr_star(:,:)=0.
     612      alim_star(:,:)=0.
     613      alim_star_tot(:)=0.
     614      csc(:,:)=0.
     615      detr(:,:)=0.
     616      entr(:,:)=0.
     617      zw2(:,:)=0.
     618      zbuoy(:,:)=0.
     619      gamma(:,:)=0.
     620      zeps(:,:)=0.
     621      w_est(:,:)=0.
     622      f_star(:,:)=0.
     623      wa_moy(:,:)=0.
     624      linter(:)=1.
     625!     linter(:)=1.
     626      b1=2.
     627! Initialisation des variables entieres
     628      lmix(:)=1
     629      lmix_bis(:)=2
     630      wmaxa(:)=0.
     631      lalim(:)=1
     632
     633!      print*,'THERMCELL PLUME ARNAUD DEDANS'
     634
     635!-------------------------------------------------------------------------
     636! On ne considere comme actif que les colonnes dont les deux premieres
     637! couches sont instables.
     638!-------------------------------------------------------------------------
     639      active(:)=ztv(:,1)>ztv(:,2)
     640
     641!-------------------------------------------------------------------------
     642! Definition de l'alimentation a l'origine dans thermcell_init
     643!-------------------------------------------------------------------------
     644      do l=1,klev-1
     645         do ig=1,ngrid
     646            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
     647               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
     648     &                       *sqrt(zlev(ig,l+1))
     649               lalim(:)=l+1
     650               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     651            endif
     652         enddo
     653      enddo
     654      do l=1,klev
     655         do ig=1,ngrid
     656            if (alim_star_tot(ig) > 1.e-10 ) then
     657               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
     658            endif
     659         enddo
     660      enddo
     661      alim_star_tot(:)=1.
     662
     663
     664
     665!------------------------------------------------------------------------------
     666! Calcul dans la premiere couche
     667! On decide dans cette version que le thermique n'est actif que si la premiere
     668! couche est instable.
     669! Pourrait etre change si on veut que le thermiques puisse se déclencher
     670! dans une couche l>1
     671!------------------------------------------------------------------------------
     672do ig=1,ngrid
     673! Le panache va prendre au debut les caracteristiques de l'air contenu
     674! dans cette couche.
     675    if (active(ig)) then
     676    ztla(ig,1)=zthl(ig,1)
     677    zqta(ig,1)=po(ig,1)
     678    zqla(ig,1)=zl(ig,1)
     679!cr: attention, prise en compte de f*(1)=1
     680    f_star(ig,2)=alim_star(ig,1)
     681    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
     682&                     *(zlev(ig,2)-zlev(ig,1))  &
     683&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
     684    w_est(ig,2)=zw2(ig,2)
     685    endif
     686enddo
     687!
     688
     689!==============================================================================
     690!boucle de calcul de la vitesse verticale dans le thermique
     691!==============================================================================
     692do l=2,klev-1
     693!==============================================================================
     694
     695
     696! On decide si le thermique est encore actif ou non
     697! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     698    do ig=1,ngrid
     699       active(ig)=active(ig) &
     700&                 .and. zw2(ig,l)>1.e-10 &
     701&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     702    enddo
     703
     704
     705
     706!---------------------------------------------------------------------------
     707! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     708! sans tenir compte du detrainement et de l'entrainement dans cette
     709! couche
     710! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     711! avant) a l'alimentation pour avoir un calcul plus propre
     712!---------------------------------------------------------------------------
     713
     714     call thermcell_condens(ngrid,active, &
     715&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
     716
     717
     718    do ig=1,ngrid
     719!       print*,'active',active(ig),ig,l
     720        if(active(ig)) then
     721        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
     722        zta_est(ig,l)=ztva_est(ig,l)
     723        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
     724        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
     725     &      -zqla_est(ig,l))-zqla_est(ig,l))
     726
     727!------------------------------------------------
     728!AJAM:nouveau calcul de w² 
     729!------------------------------------------------
     730              zdz=zlev(ig,l+1)-zlev(ig,l)
     731              zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     732
     733              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     734              zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
     735              w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
     736 
     737
     738             if (w_est(ig,l+1).lt.0.) then
     739                w_est(ig,l+1)=zw2(ig,l)
     740             endif
     741       endif
     742    enddo
     743
     744
     745!-------------------------------------------------
     746!calcul des taux d'entrainement et de detrainement
     747!-------------------------------------------------
     748
     749     do ig=1,ngrid
     750        if (active(ig)) then
     751
     752          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
     753          zw2m=w_est(ig,l+1)
     754          zdz=zlev(ig,l+1)-zlev(ig,l)
     755          zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     756!          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
     757          zbuoybis=zbuoy(ig,l)
     758          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
     759          zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
     760
     761         
     762          entr_star(ig,l)=f_star(ig,l)*zdz*  zbetalpha*MAX(0.,  &
     763    &     afact*zbuoybis/zw2m - fact_epsilon )
     764
     765
     766          detr_star(ig,l)=f_star(ig,l)*zdz                        &
     767    &     *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m          &
     768    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
     769         
     770! En dessous de lalim, on prend le max de alim_star et entr_star pour
     771! alim_star et 0 sinon
     772        if (l.lt.lalim(ig)) then
     773          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
     774          entr_star(ig,l)=0.
     775        endif
     776
     777! Calcul du flux montant normalise
     778      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
     779     &              -detr_star(ig,l)
     780
     781      endif
     782   enddo
     783
     784
     785!----------------------------------------------------------------------------
     786!calcul de la vitesse verticale en melangeant Tl et qt du thermique
     787!---------------------------------------------------------------------------
     788   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
     789   do ig=1,ngrid
     790       if (activetmp(ig)) then
     791           Zsat=.false.
     792           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
     793     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
     794     &            /(f_star(ig,l+1)+detr_star(ig,l))
     795           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
     796     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
     797     &            /(f_star(ig,l+1)+detr_star(ig,l))
     798
     799        endif
     800    enddo
     801
     802   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
     803
     804
     805   do ig=1,ngrid
     806      if (activetmp(ig)) then
     807        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
     808! on ecrit de maniere conservative (sat ou non)
     809!          T = Tl +Lv/Cp ql
     810           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
     811           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
     812!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
     813           zha(ig,l) = ztva(ig,l)
     814           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
     815     &              -zqla(ig,l))-zqla(ig,l))
     816
     817!on ecrit zqsat
     818           zqsatth(ig,l)=qsatbef 
     819
     820           zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     821           zdz=zlev(ig,l+1)-zlev(ig,l)
     822           zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
     823
     824            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
     825            zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
     826            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
     827      endif
     828   enddo
     829
     830        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
     831!
     832!---------------------------------------------------------------------------
     833!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
     834!---------------------------------------------------------------------------
     835
     836   do ig=1,ngrid
     837            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
     838!               stop'On tombe sur le cas particulier de thermcell_dry'
     839                print*,'On tombe sur le cas particulier de thermcell_plume'
     840                zw2(ig,l+1)=0.
     841                linter(ig)=l+1
     842            endif
     843
     844        if (zw2(ig,l+1).lt.0.) then
     845           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
     846     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
     847           zw2(ig,l+1)=0.
     848        endif
     849
     850           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
     851
     852        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
     853!   lmix est le niveau de la couche ou w (wa_moy) est maximum
     854!on rajoute le calcul de lmix_bis
     855            if (zqla(ig,l).lt.1.e-10) then
     856               lmix_bis(ig)=l+1
     857            endif
     858            lmix(ig)=l+1
     859            wmaxa(ig)=wa_moy(ig,l+1)
     860        endif
     861   enddo
     862
     863!=========================================================================
     864! FIN DE LA BOUCLE VERTICALE
     865      enddo
     866!=========================================================================
     867
     868!      print*,'THERMCELL PLUME ARNAUD DEDANS 7'
     869!on recalcule alim_star_tot
     870       do ig=1,ngrid
     871          alim_star_tot(ig)=0.
     872       enddo
     873       do ig=1,ngrid
     874          do l=1,lalim(ig)-1
     875          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     876          enddo
     877       enddo
     878       
     879
     880        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
     881
     882#undef wrgrads_thermcell
     883#ifdef wrgrads_thermcell
     884         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
     885         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
     886         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
     887         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
     888         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
     889         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
     890         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
     891#endif
     892
     893
     894!      print*,'THERMCELL PLUME ARNAUD DEDANS 8'
     895     return
     896     end
     897
Note: See TracChangeset for help on using the changeset viewer.