! ! SUBROUTINE thermcell_main_mars(ngrid,nlay,nq,ptimestep & & ,pplay,pplev,pphi,zlev,zlay & & ,pu,pv,pt,pq,pq2 & & ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj & & ,fm,entr,detr,lmax & & ,r_aspect & & ,zw2,fraca & & ,zpopsk,ztla,heatFlux,heatFlux_down & & ,buoyancyOut, buoyancyEst) USE thermcell,ONLY:RG,RD,RKAPPA IMPLICIT NONE !======================================================================= ! Mars version of thermcell_main. Author : A Colaitis !======================================================================= ! ============== INPUTS ============== INTEGER, INTENT(IN) :: ngrid,nlay,nq REAL, INTENT(IN) :: ptimestep,r_aspect REAL, INTENT(IN) :: pt(ngrid,nlay) REAL, INTENT(IN) :: pu(ngrid,nlay) REAL, INTENT(IN) :: pv(ngrid,nlay) REAL, INTENT(IN) :: pq(ngrid,nlay,nq) REAL, INTENT(IN) :: pq2(ngrid,nlay) REAL, INTENT(IN) :: pplay(ngrid,nlay) REAL, INTENT(IN) :: pplev(ngrid,nlay+1) REAL, INTENT(IN) :: pphi(ngrid,nlay) REAL, INTENT(IN) :: zlay(ngrid,nlay) REAL, INTENT(IN) :: zlev(ngrid,nlay+1) ! ============== OUTPUTS ============== REAL, INTENT(OUT) :: pdtadj(ngrid,nlay) REAL, INTENT(OUT) :: pduadj(ngrid,nlay) REAL, INTENT(OUT) :: pdvadj(ngrid,nlay) REAL, INTENT(OUT) :: pdqadj(ngrid,nlay,nq) REAL, INTENT(OUT) :: pdq2adj(ngrid,nlay) REAL, INTENT(OUT) :: zw2(ngrid,nlay+1) ! Diagnostics REAL, INTENT(OUT) :: heatFlux(ngrid,nlay) ! interface heatflux REAL, INTENT(OUT) :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft REAL, INTENT(OUT) :: buoyancyOut(ngrid,nlay) ! interlayer buoyancy term REAL, INTENT(OUT) :: buoyancyEst(ngrid,nlay) ! interlayer estimated buoyancy term ! dummy variables when output not needed : ! REAL :: heatFlux(ngrid,nlay) ! interface heatflux ! REAL :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft ! REAL :: buoyancyOut(ngrid,nlay) ! interlayer buoyancy term ! REAL :: buoyancyEst(ngrid,nlay) ! interlayer estimated buoyancy term ! ============== LOCAL ================ INTEGER ig,k,l,ll,iq INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid) INTEGER lmix(ngrid) INTEGER lmix_bis(ngrid) REAL linter(ngrid) REAL zmix(ngrid) REAL zmax(ngrid) REAL ztva(ngrid,nlay),zw_est(ngrid,nlay+1),ztva_est(ngrid,nlay) REAL zmax_sec(ngrid) REAL zh(ngrid,nlay) REAL zdthladj(ngrid,nlay) REAL zdthladj_down(ngrid,nlay) REAL ztvd(ngrid,nlay) REAL ztv(ngrid,nlay) REAL zu(ngrid,nlay),zv(ngrid,nlay),zo(ngrid,nlay) REAL zva(ngrid,nlay) REAL zua(ngrid,nlay) REAL zta(ngrid,nlay) REAL fraca(ngrid,nlay+1) REAL q2(ngrid,nlay) REAL rho(ngrid,nlay),rhobarz(ngrid,nlay),masse(ngrid,nlay) REAL zpopsk(ngrid,nlay) REAL wmax(ngrid) REAL wmax_sec(ngrid) REAL fm(ngrid,nlay+1),entr(ngrid,nlay),detr(ngrid,nlay) REAL fm_down(ngrid,nlay+1) REAL ztla(ngrid,nlay) REAL f_star(ngrid,nlay+1),entr_star(ngrid,nlay) REAL detr_star(ngrid,nlay) REAL alim_star_tot(ngrid) REAL alim_star(ngrid,nlay) REAL alim_star_clos(ngrid,nlay) REAL f(ngrid) REAL teta_th_int(ngrid,nlay) REAL teta_env_int(ngrid,nlay) REAL teta_down_int(ngrid,nlay) CHARACTER (LEN=20) :: modname CHARACTER (LEN=80) :: abort_message ! ============= PLUME VARIABLES ============ REAL w_est(ngrid,nlay+1) REAL wa_moy(ngrid,nlay+1) REAL wmaxa(ngrid) REAL zdz,zbuoy(ngrid,nlay),zw2m LOGICAL active(ngrid),activetmp(ngrid) ! ========================================== ! ============= HEIGHT VARIABLES =========== REAL num(ngrid) REAL denom(ngrid) REAL zlevinter(ngrid) ! ========================================= ! ============= DRY VARIABLES ============= REAL zw2_dry(ngrid,nlay+1) REAL f_star_dry(ngrid,nlay+1) REAL ztva_dry(ngrid,nlay+1) REAL wmaxa_dry(ngrid) REAL wa_moy_dry(ngrid,nlay+1) REAL linter_dry(ngrid),zlevinter_dry(ngrid) INTEGER lmix_dry(ngrid),lmax_dry(ngrid) ! ========================================= ! ============= CLOSURE VARIABLES ========= REAL zdenom(ngrid) REAL alim_star2(ngrid) REAL alim_star_tot_clos(ngrid) INTEGER llmax ! ========================================= ! ============= FLUX2 VARIABLES =========== INTEGER ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha INTEGER ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8 REAL zfm REAL f_old,ddd0,eee0,ddd,eee,zzz REAL fomass_max,alphamax ! ========================================= ! ============= DTETA VARIABLES =========== ! rien : on prends la divergence du flux turbulent ! ========================================= ! ============= DV2 VARIABLES ============= ! not used for now REAL qa(ngrid,nlay),detr_dv2(ngrid,nlay),zf,zf2 REAL wvd(ngrid,nlay+1),wud(ngrid,nlay+1) REAL gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1) REAL ue(ngrid,nlay),ve(ngrid,nlay) LOGICAL ltherm(ngrid,nlay) REAL dua(ngrid,nlay),dva(ngrid,nlay) INTEGER iter INTEGER nlarga0 ! ========================================= !----------------------------------------------------------------------- ! initialisation: ! --------------- zu(:,:)=pu(:,:) zv(:,:)=pv(:,:) ztv(:,:)=pt(:,:)/zpopsk(:,:) !------------------------------------------------------------------------ ! -------------------- ! ! ! + + + + + + + + + + + ! ! ! 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(:,nlay)-pphi(:,nlay-1))/RG ! zlay(:,:)=pphi(:,:)/RG !----------------------------------------------------------------------- ! Calcul des densites !----------------------------------------------------------------------- rho(:,:)=pplay(:,:)/(RD*pt(:,:)) rhobarz(:,1)=rho(:,1) 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 !------------------------------------------------------------------ ! ! /|\ ! -------- | 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 a la base du ! panache thermique, calcule a partir de la flotabilite 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. lmin=1 !----------------------------------------------------------------------------- ! 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 utilise pour determiner la geometrie du thermique. !------------------------------------------------------------------------------ !--------------------------------------------------------------------------------- !calcul du melange et des variables dans le thermique !-------------------------------------------------------------------------------- ! =========================================================================== ! ===================== PLUME =============================================== ! =========================================================================== ! Initialisations des variables reeles ztva(:,:)=ztv(:,:) ztva_est(:,:)=ztva(:,:) ztla(:,:)=0. zdz=0. zbuoy(:,:)=0. w_est(:,:)=0. f_star(:,:)=0. wa_moy(:,:)=0. linter(:)=1. ! Initialisation des variables entieres lmix(:)=1 lmix_bis(:)=2 wmaxa(:)=0. lalim(:)=1 !------------------------------------------------------------------------- ! On ne considere comme actif que les colonnes dont les deux premieres ! couches sont instables. !------------------------------------------------------------------------- active(:)=ztv(:,1)>ztv(:,2) do ig=1,ngrid if (ztv(ig,1)>=(ztv(ig,2))) then alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.) & & *sqrt(zlev(ig,2)) ! & *zlev(ig,2) lalim(ig)=2 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1) endif enddo do l=2,nlay-1 do ig=1,ngrid if (ztv(ig,l)>(ztv(ig,l+1)+0.5) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & & *sqrt(zlev(ig,l+1)) ! & *zlev(ig,2) lalim(ig)=l+1 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) endif enddo enddo do l=1,nlay do ig=1,ngrid if (alim_star_tot(ig) > 1.e-10 ) then alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) endif enddo enddo alim_star_tot(:)=1. if(alim_star(1,1) .ne. 0.) then print*, alim_star(:,:) endif !------------------------------------------------------------------------------ ! Calcul dans la premiere couche ! On decide dans cette version que le thermique n'est actif que si la premiere ! couche est instable. ! Pourrait etre change si on veut que le thermiques puisse se déclencher ! dans une couche l>1 !------------------------------------------------------------------------------ do ig=1,ngrid ! Le panache va prendre au debut les caracteristiques de l'air contenu ! dans cette couche. if (active(ig)) then ztla(ig,1)=ztv(ig,1) !cr: attention, prise en compte de f*(1)=1 => AC : what ? f*(1) =0. ! (d'ou f*(2)=a*(1) ! dans un panache conservatif f_star(ig,1)=0. f_star(ig,2)=alim_star(ig,1) zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & & *(zlev(ig,2)-zlev(ig,1)) & & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) w_est(ig,2)=zw2(ig,2) endif enddo !============================================================================== !boucle de calcul de la vitesse verticale dans le thermique !============================================================================== do l=2,nlay-1 !============================================================================== ! On decide si le thermique est encore actif ou non ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test do ig=1,ngrid active(ig)=active(ig) & & .and. zw2(ig,l)>1.e-10 & & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 enddo !--------------------------------------------------------------------------- ! calcul des proprietes thermodynamiques et de la vitesse de la couche l ! sans tenir compte du detrainement et de l'entrainement dans cette ! couche ! C'est a dire qu'on suppose ! ztla(l)=ztla(l-1) ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer ! avant) a l'alimentation pour avoir un calcul plus propre !--------------------------------------------------------------------------- do ig=1,ngrid if(active(ig)) then ! if(l .lt. lalim(ig)) then ! ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & ! & alim_star(ig,l)*ztv(ig,l)) & ! & /(f_star(ig,l)+alim_star(ig,l)) ! else ztva_est(ig,l)=ztla(ig,l-1) ! endif zdz=zlev(ig,l+1)-zlev(ig,l) zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) if (((2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015 & & -2.*zdz*w_est(ig,l)*0.045*(2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015)**0.6) else w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015) endif if (w_est(ig,l+1).lt.0.) then w_est(ig,l+1)=zw2(ig,l) endif endif enddo !------------------------------------------------- !calcul des taux d'entrainement et de detrainement !------------------------------------------------- do ig=1,ngrid if (active(ig)) then zw2m=w_est(ig,l+1) if((2.5*(zbuoy(ig,l)/zw2m)-0.0015) .gt. 0.) then entr_star(ig,l)=f_star(ig,l)*zdz* & ! & 0.0118*((zbuoy(ig,l)/zw2m)/0.043)**(1./1.65) ! & 0.0090*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.90) ! & 0.0120*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.6) & MAX(0.,0.045*(2.5*(zbuoy(ig,l)/zw2m)-0.0015)**0.6) else entr_star(ig,l)=0. endif if(zbuoy(ig,l) .gt. 0.) then if(l .lt. lalim(ig)) then detr_star(ig,l)=0. else ! detr_star(ig,l) = f_star(ig,l)*zdz* & ! & 0.0105*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.7) ! detr_star(ig,l) = f_star(ig,l)*zdz* & ! & 0.0085*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.55) ! last baseline from direct les ! detr_star(ig,l) = f_star(ig,l)*zdz* & ! & 0.065*(2.5*(zbuoy(ig,l)/zw2m))**0.75 ! new param from continuity eq with a fit on dfdz detr_star(ig,l) = f_star(ig,l)*zdz* & & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) ! & 0.014*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.35) ! detr_star(ig,l) = f_star(ig,l)*zdz* & ! & ((zbuoy(ig,l)/zw2m)/2.222)! + 0.0002) endif else detr_star(ig,l)=f_star(ig,l)*zdz* & & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) ! & *5.*(-afact*zbetalpha*zbuoy(ig,l)/zw2m) ! & *5.*(-afact*zbuoy(ig,l)/zw2m) ! last baseline from direct les ! & 0.065*(-2.5*(zbuoy(ig,l)/zw2m))**0.75 ! new param from continuity eq with a fit on dfdz endif ! En dessous de lalim, on prend le max de alim_star et entr_star pour ! alim_star et 0 sinon if (l.lt.lalim(ig)) then alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) entr_star(ig,l)=0. endif ! Calcul du flux montant normalise f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & & -detr_star(ig,l) endif enddo !---------------------------------------------------------------------------- !calcul de la vitesse verticale en melangeant Tl et qt du thermique !--------------------------------------------------------------------------- activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 do ig=1,ngrid if (activetmp(ig)) then ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & & (alim_star(ig,l)+entr_star(ig,l))*ztv(ig,l)) & & /(f_star(ig,l+1)+detr_star(ig,l)) endif enddo do ig=1,ngrid if (activetmp(ig)) then ztva(ig,l) = ztla(ig,l) zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) if (((2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*2.5*zbuoy(ig,l)- & & 2.*zdz*zw2(ig,l)*0.0015-2.*zdz*zw2(ig,l)*0.045*(2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015)**0.6) else zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*0.0015) endif endif enddo !--------------------------------------------------------------------------- !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max !--------------------------------------------------------------------------- do ig=1,ngrid if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then print*,'On tombe sur le cas particulier de thermcell_plume' zw2(ig,l+1)=0. linter(ig)=l+1 endif if (zw2(ig,l+1).lt.0.) then linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) zw2(ig,l+1)=0. endif wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) if (wa_moy(ig,l+1).gt.wmaxa(ig)) then ! lmix est le niveau de la couche ou w (wa_moy) est maximum !on rajoute le calcul de lmix_bis lmix(ig)=l+1 wmaxa(ig)=wa_moy(ig,l+1) endif enddo !========================================================================= ! FIN DE LA BOUCLE VERTICALE enddo !========================================================================= !on recalcule alim_star_tot do ig=1,ngrid alim_star_tot(ig)=0. enddo do ig=1,ngrid do l=1,lalim(ig)-1 alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) enddo enddo do l=1,nlay do ig=1,ngrid if (alim_star_tot(ig) > 1.e-10 ) then alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) endif enddo enddo ! =========================================================================== ! ================= FIN PLUME =============================================== ! =========================================================================== ! =========================================================================== ! ================= HEIGHT ================================================== ! =========================================================================== ! Attention, w2 est transforme en sa racine carree dans cette routine !------------------------------------------------------------------------------- ! Calcul des caracteristiques du thermique:zmax,zmix,wmax !------------------------------------------------------------------------------- !calcul de la hauteur max du thermique do ig=1,ngrid lmax(ig)=lalim(ig) enddo do ig=1,ngrid do l=nlay,lalim(ig)+1,-1 if (zw2(ig,l).le.1.e-10) then lmax(ig)=l-1 endif enddo enddo ! On traite le cas particulier qu'il faudrait éviter ou le thermique ! atteind le haut du modele ... do ig=1,ngrid if ( zw2(ig,nlay) > 1.e-10 ) then print*,'WARNING !!!!! W2 thermiques non nul derniere couche ' lmax(ig)=nlay endif enddo ! pas de thermique si couche 1 stable ! do ig=1,ngrid ! if (lmin(ig).gt.1) then ! lmax(ig)=1 ! lmin(ig)=1 ! lalim(ig)=1 ! endif ! enddo ! ! Determination de zw2 max do ig=1,ngrid wmax(ig)=0. enddo do l=1,nlay do ig=1,ngrid if (l.le.lmax(ig)) then if (zw2(ig,l).lt.0.)then print*,'pb2 zw2<0' endif zw2(ig,l)=sqrt(zw2(ig,l)) wmax(ig)=max(wmax(ig),zw2(ig,l)) else zw2(ig,l)=0. endif enddo enddo ! Longueur caracteristique correspondant a la hauteur des thermiques. do ig=1,ngrid zmax(ig)=0. zlevinter(ig)=zlev(ig,1) enddo num(:)=0. denom(:)=0. do ig=1,ngrid do l=1,nlay num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) enddo enddo do ig=1,ngrid if (denom(ig).gt.1.e-10) then zmax(ig)=2.*num(ig)/denom(ig) endif enddo ! def de zmix continu (profil parabolique des vitesses) do ig=1,ngrid if (lmix(ig).gt.1) then ! test if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & & *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) & & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10) & & then ! zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & & *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2) & & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & & *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2)) & & /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & & *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) & & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) else zmix(ig)=zlev(ig,lmix(ig)) print*,'pb zmix' endif else zmix(ig)=0. endif !test if ((zmax(ig)-zmix(ig)).le.0.) then zmix(ig)=0.9*zmax(ig) endif enddo ! ! calcul du nouveau lmix correspondant do ig=1,ngrid do l=1,nlay if (zmix(ig).ge.zlev(ig,l).and. & & zmix(ig).lt.zlev(ig,l+1)) then lmix(ig)=l endif enddo enddo ! Attention, w2 est transforme en sa racine carree dans cette routine ! =========================================================================== ! ================= FIN HEIGHT ============================================== ! =========================================================================== ! =========================================================================== ! ================= DRY ===================================================== ! =========================================================================== !initialisations do ig=1,ngrid do l=1,nlay+1 zw2_dry(ig,l)=0. wa_moy_dry(ig,l)=0. enddo enddo do ig=1,ngrid do l=1,nlay ztva_dry(ig,l)=ztv(ig,l) enddo enddo do ig=1,ngrid wmax_sec(ig)=0. wmaxa_dry(ig)=0. enddo !calcul de la vitesse a partir de la CAPE en melangeant thetav ! Calcul des F^*, integrale verticale de E^* f_star_dry(:,1)=0. do l=1,nlay f_star_dry(:,l+1)=f_star_dry(:,l)+alim_star(:,l) enddo ! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise linter_dry(:)=0. ! couche la plus haute concernee par le thermique. lmax_dry(:)=1 ! Le niveau linter est une variable continue qui se trouve dans la couche ! lmax do l=1,nlay-2 do ig=1,ngrid if (l.eq.lmin(ig).and.lalim(ig).gt.1) then !------------------------------------------------------------------------ ! Calcul de la vitesse en haut de la premiere couche instable. ! Premiere couche du panache thermique !------------------------------------------------------------------------ zw2_dry(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) & & *(zlev(ig,l+1)-zlev(ig,l)) & & *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l)) !------------------------------------------------------------------------ ! Tant que la vitesse en bas de la couche et la somme du flux de masse ! et de l'entrainement (c'est a dire le flux de masse en haut) sont ! positifs, on calcul ! 1. le flux de masse en haut f_star(ig,l+1) ! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l) ! 3. la vitesse au carré en haut zw2(ig,l+1) !------------------------------------------------------------------------ else if (zw2_dry(ig,l).ge.1e-10) then ztva_dry(ig,l)=(f_star_dry(ig,l)*ztva_dry(ig,l-1)+alim_star(ig,l) & & *ztv(ig,l))/f_star_dry(ig,l+1) zw2_dry(ig,l+1)=zw2_dry(ig,l)*(f_star_dry(ig,l)/f_star_dry(ig,l+1))**2+ & & 2.*RG*(ztva_dry(ig,l)-ztv(ig,l))/ztv(ig,l) & & *(zlev(ig,l+1)-zlev(ig,l)) endif ! determination de zmax continu par interpolation lineaire !------------------------------------------------------------------------ if (zw2_dry(ig,l+1)>0. .and. zw2_dry(ig,l+1).lt.1.e-10) then ! stop'On tombe sur le cas particulier de thermcell_dry' ! print*,'On tombe sur le cas particulier de thermcell_dry' zw2_dry(ig,l+1)=0. linter_dry(ig)=l+1 lmax_dry(ig)=l endif if (zw2_dry(ig,l+1).lt.0.) then linter_dry(ig)=(l*(zw2_dry(ig,l+1)-zw2_dry(ig,l)) & & -zw2_dry(ig,l))/(zw2_dry(ig,l+1)-zw2_dry(ig,l)) zw2_dry(ig,l+1)=0. lmax_dry(ig)=l endif wa_moy_dry(ig,l+1)=sqrt(zw2_dry(ig,l+1)) if (wa_moy_dry(ig,l+1).gt.wmaxa_dry(ig)) then ! lmix est le niveau de la couche ou w (wa_moy) est maximum lmix_dry(ig)=l+1 wmaxa_dry(ig)=wa_moy_dry(ig,l+1) endif enddo enddo ! ! Determination de zw2 max do ig=1,ngrid wmax_sec(ig)=0. enddo do l=1,nlay do ig=1,ngrid if (l.le.lmax_dry(ig)) then zw2_dry(ig,l)=sqrt(zw2_dry(ig,l)) wmax_sec(ig)=max(wmax_sec(ig),zw2_dry(ig,l)) else zw2_dry(ig,l)=0. endif enddo enddo ! Longueur caracteristique correspondant a la hauteur des thermiques. do ig=1,ngrid zmax_sec(ig)=0. zlevinter_dry(ig)=zlev(ig,1) enddo do ig=1,ngrid ! calcul de zlevinter zlevinter_dry(ig)=zlev(ig,lmax_dry(ig)) + & & (linter_dry(ig)-lmax_dry(ig))*(zlev(ig,lmax_dry(ig)+1)-zlev(ig,lmax_dry(ig))) zmax_sec(ig)=max(zmax_sec(ig),zlevinter_dry(ig)-zlev(ig,lmin(ig))) enddo ! =========================================================================== ! ============= FIN DRY ===================================================== ! =========================================================================== ! Choix de la fonction d'alimentation utilisee pour la fermeture. alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:) ! =========================================================================== ! ============= CLOSURE ===================================================== ! =========================================================================== !------------------------------------------------------------------------------- ! Fermeture,determination de f !------------------------------------------------------------------------------- ! Appel avec la version seche alim_star2(:)=0. alim_star_tot_clos(:)=0. f(:)=0. ! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine llmax=1 do ig=1,ngrid if (lalim(ig)>llmax) llmax=lalim(ig) enddo ! Calcul des integrales sur la verticale de alim_star et de ! alim_star^2/(rho dz) do k=1,llmax-1 do ig=1,ngrid if (k1.e-10) then f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/ & & (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig)) endif enddo ! =========================================================================== ! ============= FIN CLOSURE ================================================= ! =========================================================================== ! =========================================================================== ! ============= FLUX2 ======================================================= ! =========================================================================== !------------------------------------------------------------------------------- !deduction des flux !------------------------------------------------------------------------------- fomass_max=0.8 alphamax=0.5 ncorecfm1=0 ncorecfm2=0 ncorecfm3=0 ncorecfm4=0 ncorecfm5=0 ncorecfm6=0 ncorecfm7=0 ncorecfm8=0 ncorecalpha=0 !------------------------------------------------------------------------- ! Multiplication par le flux de masse issu de la femreture !------------------------------------------------------------------------- do l=1,nlay entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l)) detr(:,l)=f(:)*detr_star(:,l) enddo do l=1,nlay do ig=1,ngrid if (l.lt.lmax(ig)) then fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l) elseif(l.eq.lmax(ig)) then fm(ig,l+1)=0. detr(ig,l)=fm(ig,l)+entr(ig,l) else fm(ig,l+1)=0. endif enddo enddo ! Test provisoire : pour comprendre pourquoi on corrige plein de fois ! le cas fm6, on commence par regarder une premiere fois avant les ! autres corrections. do l=1,nlay do ig=1,ngrid if (detr(ig,l).gt.fm(ig,l)) then ncorecfm8=ncorecfm8+1 endif enddo enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! FH Version en cours de test; ! par rapport a thermcell_flux, on fait une grande boucle sur "l" ! et on modifie le flux avec tous les contr�les appliques d'affilee ! pour la meme couche ! Momentanement, on duplique le calcule du flux pour pouvoir comparer ! les flux avant et apres modif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do l=1,nlay do ig=1,ngrid if (l.lt.lmax(ig)) then fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l) elseif(l.eq.lmax(ig)) then fm(ig,l+1)=0. detr(ig,l)=fm(ig,l)+entr(ig,l) else fm(ig,l+1)=0. endif enddo !------------------------------------------------------------------------- ! Verification de la positivite des flux de masse !------------------------------------------------------------------------- do ig=1,ngrid if (fm(ig,l+1).lt.0.) then print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1) ncorecfm1=ncorecfm1+1 fm(ig,l+1)=fm(ig,l) detr(ig,l)=entr(ig,l) endif enddo ! Les "optimisations" du flux sont desactivees : moins de bidouilles ! je considere que celles ci ne sont pas justifiees ou trop delicates ! pour MARS, d'apres des observations LES. !------------------------------------------------------------------------- !Test sur fraca croissant !------------------------------------------------------------------------- ! if (iflag_thermals_optflux==0) then ! do ig=1,ngrid ! if (l.ge.lalim(ig).and.l.le.lmax(ig) & ! & .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then !! zzz est le flux en l+1 a frac constant ! zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1) & ! & /(rhobarz(ig,l)*zw2(ig,l)) ! if (fm(ig,l+1).gt.zzz) then ! detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz ! fm(ig,l+1)=zzz ! ncorecfm4=ncorecfm4+1 ! endif ! endif ! enddo ! endif ! !------------------------------------------------------------------------- !test sur flux de masse croissant !------------------------------------------------------------------------- ! if (iflag_thermals_optflux==0) then ! do ig=1,ngrid ! if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then ! f_old=fm(ig,l+1) ! fm(ig,l+1)=fm(ig,l) ! detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) ! ncorecfm5=ncorecfm5+1 ! endif ! enddo ! endif ! !------------------------------------------------------------------------- !detr ne peut pas etre superieur a fm !------------------------------------------------------------------------- do ig=1,ngrid if (detr(ig,l).gt.fm(ig,l)) then ncorecfm6=ncorecfm6+1 detr(ig,l)=fm(ig,l) entr(ig,l)=fm(ig,l+1) ! Dans le cas ou on est au dessus de la couche d'alimentation et que le ! detrainement est plus fort que le flux de masse, on stope le thermique. endif if(l.gt.lmax(ig)) then detr(ig,l)=0. fm(ig,l+1)=0. entr(ig,l)=0. endif enddo !------------------------------------------------------------------------- !fm ne peut pas etre negatif !------------------------------------------------------------------------- do ig=1,ngrid if (fm(ig,l+1).lt.0.) then detr(ig,l)=detr(ig,l)+fm(ig,l+1) fm(ig,l+1)=0. ncorecfm2=ncorecfm2+1 endif enddo !----------------------------------------------------------------------- !la fraction couverte ne peut pas etre superieure a 1 !----------------------------------------------------------------------- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! FH Partie a revisiter. ! Il semble qu'etaient codees ici deux optiques dans le cas ! F/ (rho *w) > 1 ! soit limiter la hauteur du thermique en considerant que c'est ! la derniere chouche, soit limiter F a rho w. ! Dans le second cas, il faut en fait limiter a un peu moins ! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin ! dans thermcell_main et qu'il semble de toutes facons deraisonable ! d'avoir des fractions de 1.. ! Ci dessous, et dans l'etat actuel, le premier des deux if est ! sans doute inutile. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do ig=1,ngrid if (zw2(ig,l+1).gt.1.e-10) then zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax if ( fm(ig,l+1) .gt. zfm) then f_old=fm(ig,l+1) fm(ig,l+1)=zfm detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) ncorecalpha=ncorecalpha+1 endif endif enddo ! Fin de la grande boucle sur les niveaux verticaux enddo !----------------------------------------------------------------------- ! On fait en sorte que la quantite totale d'air entraine dans le ! panache ne soit pas trop grande comparee a la masse de la maille !----------------------------------------------------------------------- do l=1,nlay-1 do ig=1,ngrid eee0=entr(ig,l) ddd0=detr(ig,l) eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep ddd=detr(ig,l)-eee if (eee.gt.0.) then ncorecfm3=ncorecfm3+1 entr(ig,l)=entr(ig,l)-eee if ( ddd.gt.0.) then ! l'entrainement est trop fort mais l'exces peut etre compense par une ! diminution du detrainement) detr(ig,l)=ddd else ! l'entrainement est trop fort mais l'exces doit etre compense en partie ! par un entrainement plus fort dans la couche superieure if(l.eq.lmax(ig)) then detr(ig,l)=fm(ig,l)+entr(ig,l) else entr(ig,l+1)=entr(ig,l+1)-ddd detr(ig,l)=0. fm(ig,l+1)=fm(ig,l)+entr(ig,l) detr(ig,l)=0. endif endif endif enddo enddo ! ! ddd=detr(ig)-entre !on s assure que tout s annule bien en zmax do ig=1,ngrid fm(ig,lmax(ig)+1)=0. entr(ig,lmax(ig))=0. detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig)) enddo !----------------------------------------------------------------------- ! Impression du nombre de bidouilles qui ont ete necessaires !----------------------------------------------------------------------- !IM 090508 beg if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngrid/4. ) then print*,'thermcell warning : large number of corrections' print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',& & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', & & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', & & ncorecfm6,'x fm6', & & ncorecfm7,'x fm7', & & ncorecfm8,'x fm8', & & ncorecalpha,'x alpha' endif ! =========================================================================== ! ============= FIN FLUX2 =================================================== ! =========================================================================== ! =========================================================================== ! ============= TRANSPORT =================================================== ! =========================================================================== !------------------------------------------------------------------ ! calcul du transport vertical !------------------------------------------------------------------ ! ------------------------------------------------------------------ ! Transport de teta dans l'updraft : (remplace thermcell_dq) ! ------------------------------------------------------------------ zdthladj(:,:)=0. do ig=1,ngrid if(lmax(ig) .gt. 1) then do k=1,lmax(ig) zdthladj(ig,k)=(1./masse(ig,k))*(fm(ig,k+1)*ztv(ig,k+1)- & & fm(ig,k)*ztv(ig,k)+fm(ig,k)*ztva(ig,k)-fm(ig,k+1)*ztva(ig,k+1)) if (ztv(ig,k) + ptimestep*zdthladj(ig,k) .le. 0.) then print*,'Teta<0 in thermcell_dTeta up: qenv .. dq : ', ztv(ig,k),ptimestep*zdthladj(ig,k) if(ztv(ig,k) .gt. 0.) then zdthladj(ig,k)=0. endif endif enddo endif enddo ! ------------------------------------------------------------------ ! Prescription des proprietes du downdraft ! ------------------------------------------------------------------ ztvd(:,:)=ztv(:,:) fm_down(:,:)=0. do ig=1,ngrid if (lmax(ig) .gt. 1) then do l=1,lmax(ig) if(zlay(ig,l) .le. 0.7*zmax(ig)) then fm_down(ig,l) =-1.8*fm(ig,l) endif if(zlay(ig,l) .le. 0.06*zmax(ig)) then ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.))) elseif(zlay(ig,l) .le. 0.6*zmax(ig)) then ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.3*(ztva(ig,l)/ztv(ig,l) - 1.)) elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/0.333333)*(ztva(ig,l)/ztv(ig,l) - 1.))) else ztvd(ig,l)=ztv(ig,l) endif ! if(zlay(ig,l) .le. 0.06*zmax(ig)) then ! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig))/19.231 + 0.9938) ! elseif(zlay(ig,l) .le. 0.6*zmax(ig)) then ! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig)-0.075)/187.931 + 0.9982) ! else !! elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then ! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig) -0.60)/(-1333) + 1.00025) !! else !! ztvd(ig,l)=ztv(ig,l) ! endif enddo endif enddo ! ------------------------------------------------------------------ ! Transport de teta dans le downdraft : (remplace thermcell_dq) ! ------------------------------------------------------------------ zdthladj_down(:,:)=0. do ig=1,ngrid if(lmax(ig) .gt. 1) then do k=1,lmax(ig) zdthladj_down(ig,k)=(1./masse(ig,k))*(fm_down(ig,k+1)*ztv(ig,k+1)- & & fm_down(ig,k)*ztv(ig,k)+fm_down(ig,k)*ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1)) if (ztv(ig,k) + ptimestep*zdthladj_down(ig,k) .le. 0.) then print*,'q<0 in thermcell_dTeta down: qenv .. dq : ', ztv(ig,k),ptimestep*zdthladj_down(ig,k) if(ztv(ig,k) .gt. 0.) then zdthladj(ig,k)=0. endif endif enddo endif enddo !------------------------------------------------------------------ ! Calcul de la fraction de l'ascendance !------------------------------------------------------------------ do ig=1,ngrid fraca(ig,1)=0. fraca(ig,nlay+1)=0. enddo do l=2,nlay do ig=1,ngrid 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 ! =========================================================================== ! ============= DV2 ========================================================= ! =========================================================================== ! ROUTINE OVERIDE : ne prends pas en compte le downdraft ! de plus, le gradient de pression horizontal semble tout deregler... A VOIR if (0 .eq. 1) then !------------------------------------------------------------------ ! calcul du transport vertical du moment horizontal !------------------------------------------------------------------ ! Calcul du transport de V tenant compte d'echange par gradient ! de pression horizontal avec l'environnement ! calcul du detrainement !--------------------------- nlarga0=0. do k=1,nlay do ig=1,ngrid detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) enddo enddo ! calcul de la valeur dans les ascendances do ig=1,ngrid zua(ig,1)=zu(ig,1) zva(ig,1)=zv(ig,1) ue(ig,1)=zu(ig,1) ve(ig,1)=zv(ig,1) enddo gamma(1:ngrid,1)=0. do k=2,nlay do ig=1,ngrid ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))*ptimestep > 1.e-5*masse(ig,k) if(ltherm(ig,k).and.zmax(ig)>0.) then gamma0(ig,k)=masse(ig,k) & & *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) ) & & *0.5/zmax(ig) & & *1. else gamma0(ig,k)=0. endif if (ltherm(ig,k).and.zmax(ig)<=0.) nlarga0=nlarga0+1 enddo enddo gamma(:,:)=0. do k=2,nlay do ig=1,ngrid if (ltherm(ig,k)) then dua(ig,k)=zua(ig,k-1)-zu(ig,k-1) dva(ig,k)=zva(ig,k-1)-zv(ig,k-1) else zua(ig,k)=zu(ig,k) zva(ig,k)=zv(ig,k) ue(ig,k)=zu(ig,k) ve(ig,k)=zv(ig,k) endif enddo ! Debut des iterations !---------------------- ! AC WARNING : SALE ! do iter=1,5 do ig=1,ngrid ! Pour memoire : calcul prenant en compte la fraction reelle ! zf=0.5*(fraca(ig,k)+fraca(ig,k+1)) ! zf2=1./(1.-zf) ! Calcul avec fraction infiniement petite zf=0. zf2=1. ! la première fois on multiplie le coefficient de freinage ! par le module du vent dans la couche en dessous. ! Mais pourquoi donc ??? if (ltherm(ig,k)) then ! On choisit une relaxation lineaire. ! gamma(ig,k)=gamma0(ig,k) ! On choisit une relaxation quadratique. gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2) zua(ig,k)=(fm(ig,k)*zua(ig,k-1) & & +(zf2*entr(ig,k)+gamma(ig,k))*zu(ig,k)) & & /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2 & & +gamma(ig,k)) zva(ig,k)=(fm(ig,k)*zva(ig,k-1) & & +(zf2*entr(ig,k)+gamma(ig,k))*zv(ig,k)) & & /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2 & & +gamma(ig,k)) ! print*,' OUTPUT DV2 !!!!!!!!!!!!!!!!!!!!',k,zua(ig,k),zva(ig,k),zu(ig,k),zv(ig,k),dua(ig,k),dva(ig,k) dua(ig,k)=zua(ig,k)-zu(ig,k) dva(ig,k)=zva(ig,k)-zv(ig,k) ue(ig,k)=(zu(ig,k)-zf*zua(ig,k))*zf2 ve(ig,k)=(zv(ig,k)-zf*zva(ig,k))*zf2 endif enddo ! Fin des iterations !-------------------- enddo enddo ! k=2,nlay ! Calcul du flux vertical de moment dans l'environnement. !--------------------------------------------------------- do k=2,nlay do ig=1,ngrid wud(ig,k)=fm(ig,k)*ue(ig,k) wvd(ig,k)=fm(ig,k)*ve(ig,k) enddo enddo do ig=1,ngrid wud(ig,1)=0. wud(ig,nlay+1)=0. wvd(ig,1)=0. wvd(ig,nlay+1)=0. enddo ! calcul des tendances. !----------------------- do k=1,nlay do ig=1,ngrid pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k) & & -(entr(ig,k)+gamma(ig,k))*ue(ig,k) & & -wud(ig,k)+wud(ig,k+1)) & & /masse(ig,k) pdvadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zva(ig,k) & & -(entr(ig,k)+gamma(ig,k))*ve(ig,k) & & -wvd(ig,k)+wvd(ig,k+1)) & & /masse(ig,k) enddo enddo ! Sorties eventuelles. !---------------------- ! if (nlarga0>0) then ! print*,'WARNING !!!!!! DANS THERMCELL_DV2 ' ! print*,nlarga0,' points pour lesquels laraga=0. dans un thermique' ! print*,'Il faudrait decortiquer ces points' ! endif ! =========================================================================== ! ============= FIN DV2 ===================================================== ! =========================================================================== else modname='momentum' call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & & masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax) call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & & masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax) endif !------------------------------------------------------------------ ! incrementation dt !------------------------------------------------------------------ do l=1,nlay do ig=1,ngrid pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l) enddo enddo !------------------------------------------------------------------ ! calcul du transport vertical de traceurs !------------------------------------------------------------------ if (nq .ne. 0.) then modname='tracer' DO iq=1,nq call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & & masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax) ENDDO endif !------------------------------------------------------------------ ! calcul du transport vertical de la tke !------------------------------------------------------------------ modname='tke' call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & & masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax) ! =========================================================================== ! ============= FIN TRANSPORT =============================================== ! =========================================================================== !------------------------------------------------------------------ ! Calculs de diagnostiques pour les sorties !------------------------------------------------------------------ ! DIAGNOSTIQUE ! We compute interface values for teta env and th. The last interface ! value does not matter, as the mass flux is 0 there. do l=1,nlay-1 do ig=1,ngrid teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l)) teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l)) teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l)) enddo enddo do ig=1,ngrid teta_th_int(ig,nlay)=teta_th_int(ig,nlay-1) teta_env_int(ig,nlay)=teta_env_int(ig,nlay-1) teta_down_int(ig,nlay)=teta_down_int(ig,nlay-1) enddo do l=1,nlay do ig=1,ngrid heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l)) buoyancyOut(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) buoyancyEst(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l) enddo enddo return end