! ! $Header$ ! SUBROUTINE thermcell_main(ngrid,nlay,ptimestep & & ,pplay,pplev,pphi,debut & & ,pu,pv,pt,po & & ,pduadj,pdvadj,pdtadj,pdoadj & & ,fm0,entr0,zqla,lmax & & ,ratqscth,ratqsdiff,zqsatth & & ,r_aspect,l_mix,w2di,tho) IMPLICIT NONE !======================================================================= ! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu ! Version du 09.02.07 ! Calcul du transport vertical dans la couche limite en presence ! de "thermiques" explicitement representes avec processus nuageux ! ! Réécriture à partir d'un listing papier à Habas, le 14/02/00 ! ! le thermique est supposé homogène et dissipé par mélange avec ! son environnement. la longueur l_mix contrôle l'efficacité du ! mélange ! ! Le calcul du transport des différentes espèces se fait en prenant ! en compte: ! 1. un flux de masse montant ! 2. un flux de masse descendant ! 3. un entrainement ! 4. un detrainement ! !======================================================================= !----------------------------------------------------------------------- ! declarations: ! ------------- #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" ! arguments: ! ---------- INTEGER ngrid,nlay,w2di,tho real ptimestep,l_mix,r_aspect REAL pt(ngrid,nlay),pdtadj(ngrid,nlay) REAL pu(ngrid,nlay),pduadj(ngrid,nlay) REAL pv(ngrid,nlay),pdvadj(ngrid,nlay) REAL po(ngrid,nlay),pdoadj(ngrid,nlay) REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) real pphi(ngrid,nlay) ! local: ! ------ ! integer,save :: igout=4259 integer,save :: igout=2928 integer,save :: lunout=6 integer,save :: lev_out=10 INTEGER ig,k,l,ll real zsortie1d(klon) INTEGER lmax(klon),lmin(klon),lalim(klon) INTEGER lmix(klon) real linter(klon) real zmix(klon) real zmax(klon),zw2(klon,klev+1),ztva(klon,klev) real zmax_sec(klon) real w_est(klon,klev+1) !on garde le zmax du pas de temps precedent real zmax0(klon) save zmax0 real zlev(klon,klev+1),zlay(klon,klev) real deltaz(klon,klev) REAL zh(klon,klev),zdhadj(klon,klev) real zthl(klon,klev),zdthladj(klon,klev) REAL ztv(klon,klev) real zu(klon,klev),zv(klon,klev),zo(klon,klev) real zl(klon,klev) real zsortie(klon,klev) real zva(klon,klev) real zua(klon,klev) real zoa(klon,klev) real zta(klon,klev) real zha(klon,klev) real fraca(klon,klev+1) real zf,zf2 real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev) real q2(klon,klev) common/comtherm/thetath2,wth2 real ratqscth(klon,klev) real var real vardiff real ratqsdiff(klon,klev) integer isplit,nsplit parameter (nsplit=10) data isplit/0/ save isplit logical sorties real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev) real zpspsk(klon,klev) real wmax(klon) real wmax_sec(klon) real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev) real detr0(klon,klev) real fm(klon,klev+1),entr(klon,klev) real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev) !niveau de condensation real nivcon(klon) real zcon(klon) REAL CHI real zcon2(klon) real pcon(klon) real zqsat(klon,klev) real zqsatth(klon,klev) real f_star(klon,klev+1),entr_star(klon,klev) real detr_star(klon,klev) real alim_star_tot(klon),alim_star2(klon) real alim_star(klon,klev) real f(klon), f0(klon) save f0 real zlevinter(klon) logical debut real seuil ! character*2 str2 character*10 str10 EXTERNAL SCOPY ! !----------------------------------------------------------------------- ! initialisation: ! --------------- ! seuil=0.25 if (lev_out.ge.1) print*,'thermcell_main V4' sorties=.true. IF(ngrid.NE.klon) THEN PRINT* PRINT*,'STOP dans convadj' PRINT*,'ngrid =',ngrid PRINT*,'klon =',klon ENDIF ! !Initialisation ! do ig=1,klon if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then f0(ig)=1.e-5 zmax0(ig)=40. endif enddo !----------------------------------------------------------------------- ! Calcul de T,q,ql a partir de Tl et qT dans l environnement ! -------------------------------------------------------------------- ! CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, & & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out) if (lev_out.ge.1) print*,'thermcell_main apres thermcell_env' !------------------------------------------------------------------------ ! -------------------- ! ! ! + + + + + + + + + + + ! ! ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz ! wh,wt,wo ... ! ! + + + + + + + + + + + zh,zu,zv,zo,rho ! ! ! -------------------- zlev(1) ! \\\\\\\\\\\\\\\\\\\\ ! ! !----------------------------------------------------------------------- ! Calcul des altitudes des couches !----------------------------------------------------------------------- do l=2,nlay zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG enddo zlev(:,1)=0. zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG do l=1,nlay zlay(:,l)=pphi(:,l)/RG enddo !calcul de l epaisseur des couches do l=1,nlay deltaz(:,l)=zlev(:,l+1)-zlev(:,l) enddo ! print*,'2 OK convect8' !----------------------------------------------------------------------- ! Calcul des densites !----------------------------------------------------------------------- do l=1,nlay rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l)) enddo do l=2,nlay rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1)) enddo !calcul de la masse do l=1,nlay masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG enddo if (lev_out.ge.1) print*,'thermcell_main apres initialisation' !------------------------------------------------------------------ ! ! /|\ ! -------- | F_k+1 ------- ! ----> D_k ! /|\ <---- E_k , A_k ! -------- | F_k --------- ! ----> D_k-1 ! <---- E_k-1 , A_k-1 ! ! ! ! ! ! --------------------------- ! ! ----- F_lmax+1=0 ---------- \ ! lmax (zmax) | ! --------------------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | E ! --------------------------- | D ! | ! --------------------------- | ! | ! --------------------------- \ | ! lalim | | ! --------------------------- | | ! | | ! --------------------------- | | ! | A | ! --------------------------- | | ! | | ! --------------------------- | | ! lmin (=1 pour le moment) | | ! ----- F_lmin=0 ------------ / / ! ! --------------------------- ! ////////////////////////// ! ! !============================================================================= ! Calculs initiaux ne faisant pas intervenir les changements de phase !============================================================================= !------------------------------------------------------------------ ! 1. alim_star est le profil vertical de l'alimentation à la base du ! panache thermique, calculé à partir de la flotabilité de l'air sec ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star !------------------------------------------------------------------ ! entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0. CALL thermcell_init(ngrid,nlay,ztv,zlev, & & lalim,lmin,alim_star,alim_star_tot,lev_out) call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ') call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ') if (lev_out.ge.1) print*,'thermcell_main apres thermcell_init' if (lev_out.ge.10) then write(lunout,*) 'Dans thermcell_main 1' write(lunout,*) 'lmin ',lmin(igout) write(lunout,*) 'lalim ',lalim(igout) write(lunout,*) ' ig l alim_star thetav' write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) & & ,ztv(igout,l),l=1,lalim(igout)+4) endif !----------------------------------------------------------------------------- ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un ! panache sec conservatif (e=d=0) alimente selon alim_star ! Il s'agit d'un calcul de type CAPE ! zmax_sec est utilisé pour déterminer la géométrie du thermique. !------------------------------------------------------------------------------ ! CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & & lalim,lmin,zmax_sec,wmax_sec,lev_out) call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ') if (lev_out.ge.1) print*,'thermcell_main apres thermcell_dry' if (lev_out.ge.10) then write(lunout,*) 'Dans thermcell_main 1b' write(lunout,*) 'lmin ',lmin(igout) write(lunout,*) 'lalim ',lalim(igout) write(lunout,*) ' ig l alim_star entr_star detr_star f_star ' write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) & & ,l=1,lalim(igout)+4) endif !--------------------------------------------------------------------------------- !calcul du melange et des variables dans le thermique !-------------------------------------------------------------------------------- ! CALL thermcell_plume(ngrid,nlay,ztv,zthl,po,zl,rhobarz, & & zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star, & & lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, & & ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out) call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ') if (lev_out.ge.1) print*,'thermcell_main apres thermcell_plume' if (lev_out.ge.10) then write(lunout,*) 'Dans thermcell_main 2' write(lunout,*) 'lmin ',lmin(igout) write(lunout,*) 'lalim ',lalim(igout) write(lunout,*) ' ig l alim_star entr_star detr_star f_star ' write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) & & ,f_star(igout,l+1),l=1,nint(linter(igout))+5) endif !------------------------------------------------------------------------------- ! Calcul des caracteristiques du thermique:zmax,zmix,wmax !------------------------------------------------------------------------------- ! CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, & & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out) call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ') if (lev_out.ge.1) print*,'thermcell_main apres thermcell_height' !------------------------------------------------------------------------------- ! Fermeture,determination de f !------------------------------------------------------------------------------- CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, & & zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out) if(lev_out.ge.1)print*,'thermcell_closure apres thermcell_closure' !------------------------------------------------------------------------------- !deduction des flux !------------------------------------------------------------------------------- CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, & & lalim,lmax,alim_star, & & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, & & detr,zqla,zmax,lev_out,lunout,igout) if (lev_out.ge.1) print*,'thermcell_main apres thermcell_flux' call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ') !c------------------------------------------------------------------ ! calcul du transport vertical !------------------------------------------------------------------ if (w2di.eq.1) then fm0=fm0+ptimestep*(fm-fm0)/float(tho) entr0=entr0+ptimestep*(entr-entr0)/float(tho) else fm0=fm entr0=entr detr0=detr endif call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, & & zthl,zdthladj,zta,lev_out) call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, & & po,pdoadj,zoa,lev_out) if (1.eq.0) then ! Calcul du transport de V tenant compte d'echange par gradient ! de pression horizontal avec l'environnement call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse & & ,fraca,zmax & & ,zu,zv,pduadj,pdvadj,zua,zva,lev_out) else ! calcul purement conservatif pour le transport de V call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse & & ,zu,pduadj,zua,lev_out) call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse & & ,zv,pdvadj,zva,lev_out) endif ! print*,'13 OK convect8' do l=1,nlay do ig=1,ngrid pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) enddo enddo print*,'14 OK convect8' !------------------------------------------------------------------ ! Calculs de diagnostiques pour les sorties !------------------------------------------------------------------ !calcul de fraca pour les sorties if (sorties) then do ig=1,klon fraca(ig,1)=0. enddo do l=2,nlay do ig=1,klon if (zw2(ig,l).gt.1.e-10) then fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l)) else fraca(ig,l)=0. endif enddo enddo print*,'14a OK convect8' ! calcul du niveau de condensation ! initialisation do ig=1,ngrid nivcon(ig)=0. zcon(ig)=0. enddo !nouveau calcul do ig=1,ngrid CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1)) pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI enddo do k=1,nlay do ig=1,ngrid if ((pcon(ig).le.pplay(ig,k)) & & .and.(pcon(ig).gt.pplay(ig,k+1))) then zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100. endif enddo enddo print*,'14b OK convect8' do k=nlay,1,-1 do ig=1,ngrid if (zqla(ig,k).gt.1e-10) then nivcon(ig)=k zcon(ig)=zlev(ig,k) endif enddo enddo print*,'14c OK convect8' !calcul des moments !initialisation do l=1,nlay do ig=1,ngrid q2(ig,l)=0. wth2(ig,l)=0. wth3(ig,l)=0. ratqscth(ig,l)=0. ratqsdiff(ig,l)=0. enddo enddo print*,'14d OK convect8' do l=1,nlay do ig=1,ngrid zf=fraca(ig,l) zf2=zf/(1.-zf) thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2 wth2(ig,l)=zf2*(zw2(ig,l))**2 ! print*,'wth2=',wth2(ig,l) wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) & & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l) q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2 !test: on calcul q2/po=ratqsc ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.)) enddo enddo !calcul du ratqscdiff print*,'14e OK convect8' var=0. vardiff=0. ratqsdiff(:,:)=0. do ig=1,ngrid do l=1,lalim(ig) var=var+alim_star(ig,l)*zqta(ig,l)*1000. enddo enddo print*,'14f OK convect8' do ig=1,ngrid do l=1,lalim(ig) zf=fraca(ig,l) zf2=zf/(1.-zf) vardiff=vardiff+alim_star(ig,l) & & *(zqta(ig,l)*1000.-var)**2 ! ratqsdiff=ratqsdiff+alim_star(ig,l)* ! s (zqta(ig,l)*1000.-po(ig,l)*1000.)**2 enddo enddo print*,'14g OK convect8' do l=1,nlay do ig=1,ngrid ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.) ! write(11,*)'ratqsdiff=',ratqsdiff(ig,l) enddo enddo !-------------------------------------------------------------------- ! !ecriture des fichiers sortie ! print*,'15 OK convect8' isplit=isplit+1 #ifdef und if (lev_out.ge.1) print*,'thermcell_main sorties 1D' #include "thermcell_out1d.h" #endif ! #define troisD if (lev_out.ge.1) print*,'thermcell_main sorties 3D' #ifdef troisD #include "thermcell_out3d.h" #endif endif if (lev_out.ge.1) print*,'thermcell_main FIN OK' return end !----------------------------------------------------------------------------- subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment) integer klon,klev real pplev(klon,klev+1),pplay(klon,klev) real ztv(klon,klev) real po(klon,klev) real ztva(klon,klev) real zqla(klon,klev) real f_star(klon,klev) real zw2(klon,klev) integer long(klon) real seuil character*21 comment print*,'TEST ',comment ! test sur la hauteur des thermiques ... do i=1,klon if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then print*,'WARNING ',comment,' au point ',i,' K= ',long(i) print*,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2' do k=1,klev write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k) enddo ! stop endif enddo return end