MODULE lmdz_thermcell_down CONTAINS SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac,dtrac) USE lmdz_thermcell_ini, ONLY: iflag_thermals_down !----------------------------------------------------------------- ! thermcell_updown_dq: computes the tendency of tracers associated ! with the presence of convective up/down drafts ! This routine that has been collectively written during the ! "ateliers downdrafts" in 2022/2023 ! Maelle, Frédéric, Catherine, Fleur, Florent, Etienne !------------------------------------------------------------------ IMPLICIT NONE ! declarations !============================================================== ! input/output integer,intent(in) :: ngrid ! number of horizontal grid points integer, intent(in) :: nlay ! number of vertical layers real,intent(in) :: ptimestep ! time step of the physics [s] real,intent(in), dimension(ngrid,nlay) :: eup ! entrainment to updrafts * dz [same unit as flux] real,intent(in), dimension(ngrid,nlay) :: dup ! detrainment from updrafts * dz [same unit as flux] real,intent(in), dimension(ngrid,nlay) :: edn ! entrainment to downdrafts * dz [same unit as flux] real,intent(in), dimension(ngrid,nlay) :: ddn ! detrainment from downdrafts * dz [same unit as flux] real,intent(in), dimension(ngrid,nlay) :: masse ! mass of layers = rho dz real,intent(in), dimension(ngrid,nlay) :: trac ! tracer integer, intent(in), dimension(ngrid) :: lmax ! max level index at which downdraft are present real,intent(out),dimension(ngrid,nlay) ::dtrac ! tendance du traceur ! Local real, dimension(ngrid,nlay+1) :: fup,fdn,fc,fthu,fthd,fthe,fthtot real, dimension(ngrid,nlay) :: tracu,tracd,traci,tracold real :: www, mstar_inv integer ig,ilay real, dimension(ngrid,nlay):: s1,s2,num !coefficients pour la resolution implicite integer :: iflag_impl ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement fdn(:,:)=0. fup(:,:)=0. fc(:,:)=0. fthu(:,:)=0. fthd(:,:)=0. fthe(:,:)=0. fthtot(:,:)=0. tracd(:,:)=0. tracu(:,:)=0. traci(:,:)=trac(:,:) tracold(:,:)=trac(:,:) s1(:,:)=0. s2(:,:)=0. num(:,:)=1. iflag_impl=1 ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement if ( iflag_thermals_down < 10 ) then call abort_physic("thermcell_updown_dq", & 'thermcell_down_dq = 0 or >= 10', 1) else iflag_impl=iflag_thermals_down-10 endif ! lmax : indice tel que fu(kmax+1)=0 ! Dans ce cas, pas besoin d'initialiser tracd(lmax) ( =trac(lmax) ) ! Boucle pour le downdraft do ilay=nlay,1,-1 do ig=1,ngrid !if ( lmax(ig) > nlay - 2 ) stop "les thermiques montent trop haut" if (ilay.le.lmax(ig) .and. lmax(ig)>1 ) then fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay) if ( fdn(ig,ilay)+ddn(ig,ilay) > 0. ) then www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay)) else www=0. endif tracd(ig,ilay)=www*tracd(ig,ilay+1) + (1.-www)*trac(ig,ilay) endif enddo enddo !Fin boucle sur l'updraft fdn(:,1)=0. !Boucle pour l'updraft do ilay=1,nlay,1 do ig=1,ngrid if (ilay.lt.lmax(ig) .and. lmax(ig)>1) then fup(ig,ilay+1)=fup(ig,ilay)+eup(ig,ilay)-dup(ig,ilay) if (fup(ig,ilay+1)+dup(ig,ilay) > 0.) then www=fup(ig,ilay)/(fup(ig,ilay+1)+dup(ig,ilay)) else www=0. endif if (ilay == 1 ) then tracu(ig,ilay)=trac(ig,ilay) else tracu(ig,ilay)=www*tracu(ig,ilay-1)+(1.-www)*trac(ig,ilay) endif endif enddo enddo !fin boucle sur le downdraft ! Calcul des flux des traceurs dans les updraft et les downdrfat ! et du flux de masse compensateur ! en ilay=1 et nlay+1, fthu=0 et fthd=0 fthu(:,1)=0. fthu(:,nlay+1)=0. fthd(:,1)=0. fthd(:,nlay+1)=0. fc(:,1)=0. fc(:,nlay+1)=0. do ilay=2,nlay,1 !boucle sur les interfaces do ig=1,ngrid fthu(ig,ilay)=fup(ig,ilay)*tracu(ig,ilay-1) fthd(ig,ilay)=-fdn(ig,ilay)*tracd(ig,ilay) fc(ig,ilay)=fup(ig,ilay)-fdn(ig,ilay) enddo enddo !Boucle pour calculer le flux du traceur flux updraft, flux downdraft, flux compensatoire !Methode explicite : if(iflag_impl==0) then do ilay=2,nlay,1 do ig=1,ngrid !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!! !!!! tentative de prise en compte d'un flux compensatoire montant !!!! if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then call abort_physic("thermcell_updown_dq", 'flux compensatoire '& // 'montant, cas non traite par thermcell_updown_dq', 1) !fthe(ig,ilay)=(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1) else fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay) endif !! si on voulait le prendre en compte on !fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1) fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay) enddo enddo !Boucle pour calculer trac do ilay=1,nlay do ig=1,ngrid dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay) ! trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay)) enddo enddo !fin du calculer de la tendance du traceur avec la methode explicite !!! Reecriture du schéma explicite avec les notations du schéma implicite else if(iflag_impl==-1) then write(*,*) 'nouveau schéma explicite !!!' !!! Calcul de s1 do ilay=1,nlay do ig=1,ngrid s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1) s2(ig,ilay)=s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1) enddo enddo do ilay=2,nlay,1 do ig=1,ngrid if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then call abort_physic("thermcell_updown_dq", 'flux compensatoire ' & // 'montant, cas non traite par thermcell_updown_dq', 1) else fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay) endif fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay) enddo enddo !Boucle pour calculer trac do ilay=1,nlay do ig=1,ngrid ! dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay) dtrac(ig,ilay)=(s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))/masse(ig,ilay) ! trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay)) ! trac(ig,ilay)=trac(ig,ilay) + (s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))*(ptimestep/masse(ig,ilay)) enddo enddo !fin du calculer de la tendance du traceur avec la methode explicite else if (iflag_impl==1) then do ilay=1,nlay do ig=1,ngrid s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1) enddo enddo !Boucle pour calculer traci = trac((t+dt) do ilay=nlay-1,1,-1 do ig=1,ngrid if((fup(ig,ilay)-fdn(ig,ilay)) .lt. 0) then write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq dans le cas d une resolution implicite, ilay : ', ilay call abort_physic("thermcell_updown_dq", "", 1) else mstar_inv=ptimestep/masse(ig,ilay) traci(ig,ilay)=((traci(ig,ilay+1)*fc(ig,ilay+1)+s1(ig,ilay))*mstar_inv+tracold(ig,ilay))/(1.+fc(ig,ilay)*mstar_inv) endif enddo enddo do ilay=1,nlay do ig=1,ngrid dtrac(ig,ilay)=(traci(ig,ilay)-tracold(ig,ilay))/ptimestep enddo enddo else call abort_physic("thermcell_updown_dq", & 'valeur de iflag_impl non prevue', 1) endif RETURN END SUBROUTINE thermcell_updown_dq !========================================================================= SUBROUTINE thermcell_down(ngrid,nlay,po,pt,pu,pv,pplay,pplev, & & lmax,fup,eup,dup,theta) !-------------------------------------------------------------- !thermcell_down: calcul des propri??t??s du panache descendant. !-------------------------------------------------------------- USE lmdz_thermcell_ini, ONLY : prt_level,RLvCp,RKAPPA,RETV,fact_thermals_down IMPLICIT NONE ! arguments integer,intent(in) :: ngrid,nlay real,intent(in), dimension(ngrid,nlay) :: po,pt,pu,pv,pplay,eup,dup real,intent(in), dimension(ngrid,nlay) :: theta real,intent(in), dimension(ngrid,nlay+1) :: pplev,fup integer, intent(in), dimension(ngrid) :: lmax ! Local real, dimension(ngrid,nlay) :: edn,ddn,thetad real, dimension(ngrid,nlay+1) :: fdn integer ig,ilay real dqsat_dT logical mask(ngrid,nlay) edn(:,:)=0. ddn(:,:)=0. fdn(:,:)=0. thetad(:,:)=0. ! lmax : indice tel que fu(kmax+1)=0 ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) ) ! FH MODIFS APRES REUNIONS POUR COMMISSIONS ! quelques erreurs de declaration ! probleme si lmax=1 ce qui a l'air d'??tre le cas en d??but de simu. Devrait ??tre 0 ? ! Remarques : ! on pourrait ??crire la formule de thetad ! www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay)) ! thetad(ig,ilay)= www * thetad(ig,ilay+1) + (1.-www) * theta(ig,ilay) ! Elle a l'avantage de bien montr?? la conservation, l'id??e fondamentale dans le ! transport qu'on ne fait que sommer des "sources" au travers d'un "propagateur" ! (Green) ! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas ?? se souccier (trop) ! de la possible nulit?? du d??nominateur do ilay=nlay,1,-1 do ig=1,ngrid if (ilay.le.lmax(ig).and.lmax(ig)>1) then edn(ig,ilay)=fact_thermals_down*dup(ig,ilay) ddn(ig,ilay)=fact_thermals_down*eup(ig,ilay) fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay) thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay)) endif enddo enddo ! Suite du travail : ! Ecrire la conservervation de theta_l dans le panache descendant ! Eventuellement faire la transformation theta_l -> theta_v ! Si l'air est sec (et qu'on oublie le c??t?? theta_v) on peut ! se contenter de conserver theta. ! ! Connaissant thetadn, on peut calculer la flotabilit??. ! Connaissant la flotabilit??, on peut calculer w de proche en proche ! On peut calculer le detrainement de facon ?? garder alpha*rho = cste ! On en d??duit l'entrainement lat??ral ! C'est le mod??le des mini-projets. !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Initialisations : !------------------ ! RETURN END SUBROUTINE thermcell_down END MODULE lmdz_thermcell_down