Index: LMDZ6/trunk/libf/phylmd/thermcell_down.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/thermcell_down.F90 (revision 4436) +++ LMDZ6/trunk/libf/phylmd/thermcell_down.F90 (revision 4437) @@ -1,3 +1,6 @@ - SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac,dtrac) +SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac,dtrac) + +USE thermcell_ini_mod, ONLY: iflag_thermals_down + !----------------------------------------------------------------- @@ -32,11 +35,14 @@ ! Local - real, dimension(ngrid,nlay+1) :: fup,fdn,fthu,fthd,fthe,fthtot - real, dimension(ngrid,nlay) :: tracu,tracd - real :: www + 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=1 ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement + fdn(:,:)=0. fup(:,:)=0. + fc(:,:)=0. fthu(:,:)=0. fthd(:,:)=0. @@ -45,4 +51,16 @@ tracd(:,:)=0. tracu(:,:)=0. + traci(:,:)=trac(:,:) + tracold(:,:)=trac(:,:) + s1(:,:)=0. + s2(:,:)=0. + num(:,:)=1. + + if ( iflag_thermals_down < 10 ) then + stop 'thermcell_down_dq = 0 or >= 10' + else + iflag_impl=iflag_thermals_down-10 + endif + ! lmax : indice tel que fu(kmax+1)=0 @@ -64,5 +82,6 @@ endif enddo - enddo + enddo !Fin boucle sur l'updraft + fdn(:,1)=0. !Boucle pour l'updraft @@ -83,121 +102,115 @@ endif enddo - enddo - !Boucle pour calculer le flux du traceur flux updraft, flux downdraft, flux compensatoire - do ilay=2,nlay,1 + 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) - !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!! - !!!! si ce n'est pas le cas on stoppe le code !!!! - if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then - write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq' - stop - endif - fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay) - !! 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) + fc(ig,ilay)=fup(ig,ilay)-fdn(ig,ilay) enddo enddo - !Boucle pour calculer trac - do ilay=1,nlay,1 - do ig=1,ngrid - dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))*(1./masse(ig,ilay)) -! trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay)) - enddo - enddo -! Il reste a coder : -! d(rho trac)/dt = - d/dz(rho w'trac') -! d(trac)/dt = -1/rho * d/dz(rho w'trac') -! hydrostatique : dp - rho g dz -! dz = - dp / (rho g) -! 1 / dz = - (rho g ) / dp - -! d(trac)/dt = -1/rho * d/dp(rho w'trac') * (- rho g ) -! d(trac)/dt = d/dp(rho w'trac') * g - - -! ----> calculer d/dz(w'trac') -! w'trac' = alpha_up*(w_up-w_bar)*(trac_up-trac_bar) -! + alpha_down*(w_down-w_bar)*(trac_down-trac_bar) -! + (1-alpha_up-alpha_down)*(w_env-w_bar)*(trac_env-trac_bar) -! avec w_bar=0=alpha_up*w_up+alpha_down*w_dn+(1-alpha_up-alpha_down)w_env -! -> w_env= - (alpha_up*w_up+alpha_down*w_dn)/(1-alpha_up-alpha_down) -! et on a : trac_bar=alpha_up*trac_up + alpha_down*trac_dn + (1-alpha_up-alpha_down)trac_env -! -! w'trac' = alpha_up*w_up*(trac_up-trac_bar) -! + alpha_down*w_down*(trac_down-trac_bar) -! + (1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar) - -! rho*w'trac' = rho*alpha_up*w_up*(trac_up-trac_bar) -! + rho*alpha_down*w_down*(trac_down-trac_bar) -! + rho*(1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar) -! - -! rho*w'trac' = fup*(trac_up-trac_bar) -! + fdn*(trac_down-trac_bar) -! - rho*(alpha_up*w_up+alpha_down*w_dn)*(trac_env-trac_bar) - -! - rho*(alpha_up*w_up+alpha_down*w_dn)*((trac_bar-alpha_up*trac_up-alpha_down*trac_down) /(1-alpha_up-alpha_down)- trac_bar) -! - rho*(alpha_up*w_up+alpha_down*w_dn)*(trac_bar*(1/(1-alpha_up-alpha_down)-1) - (alpha_up*trac_up+alpha_down*trac_down) /(1-alpha_up-alpha_down)) - -! rho*w'trac' = fup*(trac_up-trac_bar) -! + fdn*(trac_down-trac_bar) -! - (fup+fdn)*(trac_env-trac_bar) - - -! rho*w'trac' = fup*(trac_up-trac_bar) -! + fdn*(trac_down-trac_bar) -! - (fup+fdn)*( (trac_bar-alpha_up*trac_up-alpha_down*trac_down) /(1-alpha_up-alpha_down) - trac_bar) - -!!! hypoth??se : alpha_up+alpha_down << 1 -> 1/(1-alpha_up-alpha_down) ~ 1 - -! rho*w'trac' = fup*(trac_up-trac_bar) - -! + fdn*(trac_down-trac_bar) -! + (fup+fdn)*(alpha_up*trac_up+alpha_down*trac_down) - -! d(trac)/dt= -1/rho d/dz(rho*w'trac') -! choix de schema temporel (euler explicite) -! (trac(t,k)-trac(t-1,k))/dt = f(t-1,k) -! (trac(t,k)-trac(t-1,k))/dt = f(t,k) ---> euler implicite -! -----> on choisit explicite en temps - -!dans le cas d'une ascendance sans downdraft : -! d/dz(rho*w'trac') = fup(k)*qa(k-1) -fup(k+1)*qa(k) - -!(rho*w'trac'(k+1)-rho*w'trac'(k))/dz = fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) - (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) -! - fdn(k) tracdn(k) - -! en continue on a : d(trac)/dt = -1/rho * d/dz(rho w'trac') -! en discretis?? on a : -! (trac(t,k)-trac(t-1,k))/dt = -1/rho * rho*w'trac'(k+1)-rho*w'trac'(k))/dz -! trac(t,k) = trac(t-1,k) - dt/rho * (rho*w'trac'(k+1)-rho*w'trac'(k))/dz -! trac(t,k) = trac(t-1,k) - dt/rho * [fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) - -! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)] - - - -! d(trac)/dt = d/dp(rho w'trac') * g -! -1/rho * d/dz(rho w'trac') = g*d/dp(rho w'trac') = g * [fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) - -! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)] - - - -! choix de sch??ma spatial -! d/dz(rho*w'trac') = (rho*w'trac'(k+1)-rho*w'trac'(k))*dz -! d/dz(rho*w'trac') = (rho*w'trac'(k)-rho*w'trac'(k-1))*dz - - - -! d/dz(rho*w'trac') = rho w'trac'(k+1)-rho w'trac'(k) - -!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -! Initialisations : -!------------------ - - -! + + + !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 + write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq' + stop + !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 + write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq' + stop + 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 + write(*,*) 'SCHEMA IMPLICITE EN COURS DE DEVELOPPEMENT !' + 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 + stop + 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 + write(*,*) 'valeur de iflag_impl non prevue' + stop + + endif + RETURN END