Changeset 1311


Ignore:
Timestamp:
Feb 18, 2010, 2:14:02 PM (15 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

Location:
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r1299 r1311  
    204204     &      ,zmax0,f0,zw2,fraca)
    205205          else if (iflag_thermals==15.or.iflag_thermals==16) then
     206
     207!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
    206208            CALL thermcell_main(itap,klon,klev,zdt  &
    207209     &      ,pplay,paprs,pphi,debut  &
    … …  
    210212     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
    211213     &      ,ratqscth,ratqsdiff,zqsatth  &
    212      &      ,r_aspect_thermals,l_mix_thermals  &
    213      &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
     214     &      ,r_aspect_thermals,l_mix_thermals &
     215     &      ,tau_thermals,iflag_thermals_ed,Ale,Alp,lalim_conv,wght_th &
    214216     &      ,zmax0,f0,zw2,fraca)
    215217           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'

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

    r1299 r1311  
    484484     
    485485       ! Initialize ok_flux_surf (for 1D model)
    486        ok_flux_surf=.FALSE.
     486       if (klon>1) ok_flux_surf=.FALSE.
    487487       
    488488       ! Initilize debug IO

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

    r1298 r1311  
    228228     $          'coherente ', i, zmasq(i), pctsrf(i, is_ter)
    229229     $          ,pctsrf(i, is_lic)
     230            zmasq(i) = fractint(i)
    230231        ENDIF
    231232      END DO
    … …  
    237238     $          'coherente ', i, zmasq(i) , pctsrf(i, is_oce)
    238239     $          ,pctsrf(i, is_sic)
     240            zmasq(i) = fractint(i)
    239241        ENDIF
    240242      END DO

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

    r1299 r1311  
    23392339con rajoute ale et alp, et les caracteristiques de la couche alim
    23402340     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
     2341
     2342! ----------------------------------------------------------------------
     2343! Transport de la TKE par les panaches thermiques.
     2344! FH : 2010/02/01
     2345      if (iflag_pbl.eq.10) then
     2346      call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
     2347     s           rg,paprs,pbl_tke)
     2348      endif
     2349! ----------------------------------------------------------------------
     2350
    23412351         endif
     2352
    23422353
    23432354
    … …  
    25502561
    25512562!   les ratqs sont une combinaison de ratqss et ratqsc
    2552        write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
     2563!       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
    25532564
    25542565         if (tau_ratqs>1.e-10) then

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

    r1279 r1311  
    8585       IF( INT( tabcntr0( 6 ) ) .NE. nbapp_rad  )   THEN
    8686        PRINT 21,  INT(tabcntr0(6)), nbapp_rad
    87         radpas0  = NINT( 86400./tabcntr0(1)/INT( tabcntr0(6) ) )
     87!        radpas0  = NINT( 86400./tabcntr0(1)/INT( tabcntr0(6) ) )
    8888        PRINT 100
    8989        PRINT 22, radpas0, radpas

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

    r1294 r1311  
    3939
    4040!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    41 print*,'THERMCELL CLOSURE 26E'
     41!print*,'THERMCELL CLOSURE 26E'
    4242
    4343alim_star2(:)=0.

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

    r1026 r1311  
    4040         enddo
    4141      enddo
     42
     43! On traite le cas particulier qu'il faudrait éviter ou le thermique
     44! atteind le haut du modele ...
     45      do ig=1,ngrid
     46      if ( zw2(ig,nlay) > 1.e-10 ) then
     47          print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
     48          lmax(ig)=nlay
     49      endif
     50      enddo
     51
    4252! pas de thermique si couche 1 stable
    4353      do ig=1,ngrid

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

    r1299 r1311  
    88     &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
    99     &                  ,ratqscth,ratqsdiff,zqsatth  &
    10      &                  ,r_aspect,l_mix,tau_thermals &
     10     &                  ,r_aspect,l_mix,tau_thermals,iflag_thermals_ed &
    1111     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
    1212     &                  ,zmax0, f0,zw2,fraca)
    … …  
    5555      INTEGER ngrid,nlay,w2di
    5656      real tau_thermals
     57      integer iflag_thermals_ed
    5758      real ptimestep,l_mix,r_aspect
    5859      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
    … …  
    228229      ENDIF
    229230!
    230      write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
     231!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
    231232     do ig=1,klon
    232233      if (prt_level.ge.20) then
    … …  
    380381      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
    381382!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    382       CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    383      &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
    384      &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
    385      &           ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    386      &            ,lev_out,lunout1,igout)
     383
     384! Gestion temporaire de plusieurs appels à thermcell_plume au travers
     385! de la variable iflag_thermals
     386
     387!      print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
     388      if (iflag_thermals_ed<=9) then
     389!         print*,'THERM NOVUELLE/NOUVELLE/ANCIENNE'
     390         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
     391     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
     392     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     393     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     394     &    ,lev_out,lunout1,igout)
     395
     396      elseif (iflag_thermals_ed<=19) then
     397! Version d'Arnaud Jam
     398!         print*,'THERM RIO et al 2010'
     399         CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
     400     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
     401     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     402     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     403     &    ,lev_out,lunout1,igout)
     404
     405      endif
     406
    387407      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
    388408
    … …  
    421441!
    422442!
    423       write(lunout,*)'THERM NOUVEAU RIO2009 31B'
     443!!      write(lunout,*)'THERM NOUVEAU XXXXX'
    424444      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
    425445    &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
    … …  
    440460
    441461
    442 print*,'THERM 26JJJ'
     462!print*,'THERM 26JJJ'
    443463
    444464! Choix de la fonction d'alimentation utilisee pour la fermeture.

  • 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

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

    r1299 r1311  
    11!
    2 ! $Id$
     2! $Header$
    33!
    44      SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
    … …  
    3838c    iflag_pbl=7 : MY 2.0.Fournier
    3939c    iflag_pbl=8 : MY 2.5
    40 c    iflag_pbl=9 : un test ?
     40c    iflag_pbl>=9 : MY 2.5 avec diffusion verticale
    4141
    4242c.......................................................................
    … …  
    6666
    6767      integer nlay,nlev
    68 cym      PARAMETER (nlay=klev)
    69 cym      PARAMETER (nlev=klev+1)
    7068
    7169      logical first
    … …  
    9896      real fl,zzz,zl0,zq2,zn2
    9997
    100 cym      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
    101 cym     s  ,lyam(klon,klev),knyam(klon,klev)
    102 cym     s  ,w2yam(klon,klev),t2yam(klon,klev)
    103       real,allocatable,save,dimension(:,:) :: rino,smyam,styam,lyam,
    104      s                                        knyam,w2yam,t2yam
    105 cym      common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam
    106 c$OMP THREADPRIVATE(rino,smyam,styam,lyam,knyam,w2yam,t2yam)
     98      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
     99     s  ,lyam(klon,klev),knyam(klon,klev)
     100     s  ,w2yam(klon,klev),t2yam(klon,klev)
    107101      logical,save :: firstcall=.true.
    108 
    109       character (len=20) :: modname='yamada4'
    110       character (len=80) :: abort_message
    111 
    112102c$OMP THREADPRIVATE(firstcall)       
    113103      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
    … …  
    123113     
    124114      if (firstcall) then
    125         allocate(rino(klon,klev+1),smyam(klon,klev),styam(klon,klev))
    126         allocate(lyam(klon,klev),knyam(klon,klev))
    127         allocate(w2yam(klon,klev),t2yam(klon,klev))
    128115        allocate(l0(klon))
    129116        firstcall=.false.
    … …  
    131118
    132119
    133       if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
    134               abort_message = 'probleme de coherence dans appel a MY'
    135               CALL abort_gcm (modname,abort_message,1)
     120      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
     121           stop'probleme de coherence dans appel a MY'
    136122      endif
    137123
    … …  
    417403      enddo
    418404
    419 c     if (iflag_pbl.ge.7..and.0.eq.1) then
    420 c        q2(:,1)=q2(:,2)
    421 c        call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
    422 c     endif
     405      if (iflag_pbl.ge.9) then
     406!       print*,'YAMADA VDIF'
     407        q2(:,1)=q2(:,2)
     408        call vdif_q2(dt,g,rconst,ngrid,plev,temp,kq,q2)
     409      endif
    423410
    424411c   Traitement des cas noctrunes avec l'introduction d'une longueur
    … …  
    497484      return
    498485      end
     486      SUBROUTINE vdif_q2(timestep,gravity,rconst,ngrid,plev,temp,
     487     &  kmy,q2)
     488      use dimphy
     489      IMPLICIT NONE
     490c.......................................................................
     491#include "dimensions.h"
     492cccc#include "dimphy.h"
     493c.......................................................................
     494c
     495c dt : pas de temps
     496
     497      real plev(klon,klev+1)
     498      real temp(klon,klev)
     499      real timestep
     500      real gravity,rconst
     501      real kstar(klon,klev+1),zz
     502      real kmy(klon,klev+1)
     503      real q2(klon,klev+1)
     504      real deltap(klon,klev+1)
     505      real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
     506      integer ngrid
     507
     508      integer i,k
     509
     510!       print*,'RD=',rconst
     511      do k=1,klev
     512         do i=1,ngrid
     513c test
     514!       print*,'i,k',i,k
     515!       print*,'temp(i,k)=',temp(i,k)
     516!       print*,'(plev(i,k)-plev(i,k+1))=',plev(i,k),plev(i,k+1)
     517            zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
     518            kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz
     519     s      /(plev(i,k)-plev(i,k+1))*timestep
     520         enddo
     521      enddo
     522
     523      do k=2,klev
     524         do i=1,ngrid
     525            deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
     526         enddo
     527      enddo
     528      do i=1,ngrid
     529         deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
     530         deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
     531         denom(i,klev+1)=deltap(i,klev+1)+kstar(i,klev)
     532         alpha(i,klev+1)=deltap(i,klev+1)*q2(i,klev+1)/denom(i,klev+1)
     533         beta(i,klev+1)=kstar(i,klev)/denom(i,klev+1)
     534      enddo
     535
     536      do k=klev,2,-1
     537         do i=1,ngrid
     538            denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))*
     539     s      kstar(i,k)+kstar(i,k-1)
     540c   correction d'un bug 10 01 2001
     541            alpha(i,k)=(q2(i,k)*deltap(i,k)
     542     s      +kstar(i,k)*alpha(i,k+1))/denom(i,k)
     543            beta(i,k)=kstar(i,k-1)/denom(i,k)
     544         enddo
     545      enddo
     546
     547c  Si on recalcule q2(1)
     548      if(1.eq.0) then
     549      do i=1,ngrid
     550         denom(i,1)=deltap(i,1)+(1-beta(i,2))*kstar(i,1)
     551         q2(i,1)=(q2(i,1)*deltap(i,1)
     552     s      +kstar(i,1)*alpha(i,2))/denom(i,1)
     553      enddo
     554      endif
     555c   sinon, on peut sauter cette boucle...
     556
     557      do k=2,klev+1
     558         do i=1,ngrid
     559            q2(i,k)=alpha(i,k)+beta(i,k)*q2(i,k-1)
     560         enddo
     561      enddo
     562
     563      return
     564      end
     565      SUBROUTINE vdif_q2e(timestep,gravity,rconst,ngrid,
     566     &   plev,temp,kmy,q2)
     567      use dimphy
     568      IMPLICIT NONE
     569c.......................................................................
     570#include "dimensions.h"
     571cccc#include "dimphy.h"
     572c.......................................................................
     573c
     574c dt : pas de temps
     575
     576      real plev(klon,klev+1)
     577      real temp(klon,klev)
     578      real timestep
     579      real gravity,rconst
     580      real kstar(klon,klev+1),zz
     581      real kmy(klon,klev+1)
     582      real q2(klon,klev+1)
     583      real deltap(klon,klev+1)
     584      real denom(klon,klev+1),alpha(klon,klev+1),beta(klon,klev+1)
     585      integer ngrid
     586
     587      integer i,k
     588
     589      do k=1,klev
     590         do i=1,ngrid
     591            zz=(plev(i,k)+plev(i,k+1))*gravity/(rconst*temp(i,k))
     592            kstar(i,k)=0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz
     593     s      /(plev(i,k)-plev(i,k+1))*timestep
     594         enddo
     595      enddo
     596
     597      do k=2,klev
     598         do i=1,ngrid
     599            deltap(i,k)=0.5*(plev(i,k-1)-plev(i,k+1))
     600         enddo
     601      enddo
     602      do i=1,ngrid
     603         deltap(i,1)=0.5*(plev(i,1)-plev(i,2))
     604         deltap(i,klev+1)=0.5*(plev(i,klev)-plev(i,klev+1))
     605      enddo
     606
     607      do k=klev,2,-1
     608         do i=1,ngrid
     609            q2(i,k)=q2(i,k)+
     610     s      ( kstar(i,k)*(q2(i,k+1)-q2(i,k))
     611     s       -kstar(i,k-1)*(q2(i,k)-q2(i,k-1)) )
     612     s      /deltap(i,k)
     613         enddo
     614      enddo
     615
     616      do i=1,ngrid
     617         q2(i,1)=q2(i,1)+
     618     s   ( kstar(i,1)*(q2(i,2)-q2(i,1))
     619     s                                      )
     620     s   /deltap(i,1)
     621         q2(i,klev+1)=q2(i,klev+1)+
     622     s   (
     623     s    -kstar(i,klev)*(q2(i,klev+1)-q2(i,klev)) )
     624     s   /deltap(i,klev+1)
     625      enddo
     626
     627      return
     628      end
Note: See TracChangeset for help on using the changeset viewer.