MODULE lmdz_thermcell_main ! $Id: lmdz_thermcell_main.F90 5512 2025-01-28 18:07:51Z ymeurdesoif $ ! ! A REGARDER !!!!!!!!!!!!!!!!! ! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023) ! ATTENTION : dans thermcell_env, on condense potentiellement de l'eau. Mais comme on ne mélange pas l'eau liquide supposant qu'il n'y en n'a pas, c'est potentiellement un souci CONTAINS subroutine thermcell_main(itap,ngrid,nlay,ptimestep & & ,pplay,pplev,pphi,debut & & ,puwind,pvwind,ptemp,p_o,ptemp_env, po_env & & ,pduadj,pdvadj,pdtadj,pdoadj & & ,fm0,entr0,detr0,zqta,zqla,lmax & & ,ratqscth,ratqsdiff,zqsatth & & ,zmax0, f0,zw2,fraca,ztv & & ,zpspsk,ztla,zthl,ztva & & ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,zcong & #ifdef ISO & ,xtpo,xtpdoadj & #endif & ) ! USE necessaires pour les lignes importees de thermcell_env USE lmdz_thermcell_ini, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level USE lmdz_thermcell_ini, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals USE lmdz_thermcell_ini, ONLY: iflag_thermals_down,fact_thermals_down USE lmdz_thermcell_ini, ONLY: iflag_thermals_tenv USE lmdz_thermcell_ini, ONLY: RD,RG USE lmdz_thermcell_down, ONLY: thermcell_updown_dq USE lmdz_thermcell_closure, ONLY: thermcell_closure USE lmdz_thermcell_dq, ONLY: thermcell_dq USE lmdz_thermcell_dry, ONLY: thermcell_dry USE lmdz_thermcell_dv2, ONLY: thermcell_dv2 USE lmdz_thermcell_env, ONLY: thermcell_env USE lmdz_thermcell_flux2, ONLY: thermcell_flux2 USE lmdz_thermcell_height, ONLY: thermcell_height USE lmdz_thermcell_plume, ONLY: thermcell_plume USE lmdz_thermcell_plume_6A, ONLY: thermcell_plume_6A,thermcell_plume_5B ! USE necessaires pour les lignes importees de thermcell_env USE lmdz_thermcell_ini, ONLY : RLvCp,RKAPPA,RETV USE lmdz_thermcell_qsat, ONLY : thermcell_qsat #ifdef ISO USE infotrac_phy, ONLY : ntiso #ifdef ISOVERIF USE isotopes_mod, ONLY : iso_eau,iso_HDO USE isotopes_verif_mod, ONLY: iso_verif_egalite, & iso_verif_aberrant_encadre #endif #endif 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 ! ! Reecriture a partir d'un listing papier a Habas, le 14/02/00 ! ! le thermique est suppose homogene et dissipe par melange avec ! son environnement. la longueur l_mix controle l'efficacite du ! melange ! ! Le calcul du transport des differentes especes se fait en prenant ! en compte: ! 1. un flux de masse montant ! 2. un flux de masse descendant ! 3. un entrainement ! 4. un detrainement ! ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr) ! Introduction of an implicit computation of vertical advection in ! the environment of thermal plumes in thermcell_dq ! impl = 0 : explicit, 1 : implicit, -1 : old version ! controled by iflag_thermals = ! 15, 16 run with impl=-1 : numerical convergence with NPv3 ! 17, 18 run with impl=1 : more stable ! 15 and 17 correspond to the activation of the stratocumulus "bidouille" ! ! Using ! abort_physic ! iso_verif_aberrant_encadre ! iso_verif_egalite ! test_ltherm ! thermcell_closure ! thermcell_dq ! thermcell_dry ! thermcell_dv2 ! thermcell_env ! thermcell_flux2 ! thermcell_height ! thermcell_plume ! thermcell_plume_5B ! thermcell_plume_6A ! !======================================================================= !----------------------------------------------------------------------- ! declarations: ! ------------- ! arguments: ! ---------- integer, intent(in) :: itap,ngrid,nlay real, intent(in) :: ptimestep real, intent(in), dimension(ngrid,nlay) :: ptemp,puwind,pvwind,pplay,pphi,ptemp_env,po_env ! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023) real, intent(in), dimension(ngrid,nlay) :: p_o real, intent(out), dimension(ngrid,nlay) :: zpspsk real, intent(in), dimension(ngrid,nlay+1) :: pplev integer, intent(out), dimension(ngrid) :: lmax real, intent(out), dimension(ngrid,nlay) :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0 real, intent(out), dimension(ngrid,nlay) :: ztla,zqla,zqta,zqsatth,zthl real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca real, intent(inout), dimension(ngrid) :: zmax0,f0 real, intent(out), dimension(ngrid,nlay) :: ztva,ztv logical, intent(in) :: debut real,intent(out), dimension(ngrid,nlay) :: ratqscth,ratqsdiff real, intent(out), dimension(ngrid) :: pcon real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3 real, intent(out), dimension(ngrid) :: wmax_sec integer,intent(out), dimension(ngrid) :: lalim real, intent(out), dimension(ngrid,nlay+1) :: fm real, intent(out), dimension(ngrid,nlay) :: alim_star real, intent(out), dimension(ngrid) :: zmax,zcong ! local: ! ------ integer, parameter :: igout=1 integer, parameter :: lunout1=6 integer, parameter :: lev_out=10 real lambda, zf,zf2,var,vardiff,CHI integer ig,k,l,ierr,ll logical sorties real, dimension(ngrid) :: linter,zmix, zmax_sec,lintercong integer,dimension(ngrid) :: lmin,lmix,lmix_bis,nivcon, lcong real, dimension(ngrid,nlay) :: ztva_est real, dimension(ngrid,nlay) :: deltaz,zlay,zdthladj,zu,zv,z_o,zl,zva,zua,z_oa real, dimension(ngrid,nlay) :: ztemp_env ! temperarure liquide de l'environnement real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2 real, dimension(ngrid,nlay) :: rho,masse real, dimension(ngrid,nlay+1) :: zw_est,zlev real, dimension(ngrid) :: wmax,wmax_tmp real, dimension(ngrid,nlay+1) :: f_star real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos real, dimension(ngrid,nlay) :: zqsat,csc real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f real, dimension(ngrid,nlay) :: entrdn,detrdn logical, dimension(ngrid,nlay) :: mask character (len=20), parameter :: modname='thermcell_main' character (len=80) :: abort_message #ifdef ISO REAL xtpo(ntiso,ngrid,nlay),xtpdoadj(ntiso,ngrid,nlay) REAL xtzo(ntiso,ngrid,nlay) REAL xtpdoadj_tmp(ngrid,nlay) REAL xtpo_tmp(ngrid,nlay) REAL xtzo_tmp(ngrid,nlay) integer ixt #endif ! !----------------------------------------------------------------------- ! initialisation: ! --------------- ! fm=0. ; entr=0. ; detr=0. if (prt_level.ge.1) print*,'thermcell_main V4' sorties=.true. IF(ngrid.NE.ngrid) THEN PRINT*,'STOP dans convadj' PRINT*,'ngrid =',ngrid PRINT*,'ngrid =',ngrid ENDIF ! !print*,'thermcell_main debut' ! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)' do ig=1,ngrid f0(ig)=max(f0(ig),1.e-2) zmax0(ig)=max(zmax0(ig),40.) !IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2 enddo if (prt_level.ge.20) then do ig=1,ngrid print*,'th_main ig f0',ig,f0(ig) enddo endif !----------------------------------------------------------------------- ! Calcul de T,q,ql a partir de Tl et qT dans l environnement ! -------------------------------------------------------------------- ! ! On condense l'eau liquide si besoin. ! En fait on arrive ici d'habitude (jusque 6A) après réévaporation ! Dans une nouvelle mouture, on passe les profiles ! avant la couche limite : iflag_thermals_tenv=1 ! dés le début de la physique : iflag_thermals_tenv=2 ! Mais même pour 2) on ne veut sans doute pas réévaporer. ! On veut comparer thetav dans le thermique, après condensation, ! avec le theta_v effectif de l'environnement. if (iflag_thermals_tenv - 10 * ( iflag_thermals_tenv / 10 ) == 0) then CALL thermcell_env(ngrid,nlay,p_o,ptemp_env,puwind,pvwind,pplay, & & pplev,z_o,ztemp_env,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lcong,lintercong,lev_out) else ! Chantier en cours : ne pas effacer (Fredho). 15 septembre 2023 ! Dans la version originale de thermcell_env, on condense l'eau de l'environnement ! pour calculer une temperature potentielle liquide. ! On en déduit un Theta v. ! ... ! contenu de thermcell_env ! SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, & ! & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out) ! contenu thermcell_env : call thermcell_qsat(ngrid, nlay,mask,pplev,pt,po,pqsat) ! contenu thermcell_env : do ll=1,nlay ! contenu thermcell_env : do ig=1,ngrid ! contenu thermcell_env : zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll)) ! contenu thermcell_env : zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll) ! T = Tl + Lv/Cp ql ! contenu thermcell_env : zo(ig,ll) = po(ig,ll)-zl(ig,ll) ! contenu thermcell_env : enddo ! contenu thermcell_env : enddo ! contenu thermcell_env : do ll=1,nlay ! contenu thermcell_env : do ig=1,ngrid ! contenu thermcell_env : zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA ! contenu thermcell_env : zu(ig,ll)=pu(ig,ll) ! contenu thermcell_env : zv(ig,ll)=pv(ig,ll) ! contenu thermcell_env : ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll) ! contenu thermcell_env : ztv(ig,ll)=ztv(ig,ll)*(1.+RETV*(zo(ig,ll))-zl(ig,ll)) ! contenu thermcell_env : zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll) ! contenu thermcell_env : enddo ! contenu thermcell_env : enddo do l=1,nlay do ig=1,ngrid zl(ig,l)=0. zu(ig,l)=puwind(ig,l) zv(ig,l)=pvwind(ig,l) ztemp_env(ig,l)=ptemp_env(ig,l) zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA ztv(ig,l)=ztemp_env(ig,l)/zpspsk(ig,l) ztv(ig,l)=ztv(ig,l)*(1.+RETV*po_env(ig,l)) zthl(ig,l)=ptemp(ig,l)/zpspsk(ig,l) mask(ig,l)=.true. enddo enddo call thermcell_qsat(ngrid, nlay, mask,pplev,ptemp_env,p_o,zqsat) endif if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env' !------------------------------------------------------------------------ ! -------------------- ! ! ! + + + + + + + + + + + ! ! ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz ! wh,wt,wo ... ! ! + + + + + + + + + + + zh,zu,zv,z_o,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 do l=1,nlay zlay(:,l)=pphi(:,l)/RG enddo do l=1,nlay deltaz(:,l)=zlev(:,l+1)-zlev(:,l) enddo !----------------------------------------------------------------------- ! Calcul des densites et masses !----------------------------------------------------------------------- rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:)) if (prt_level.ge.10) write(lunout,*) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)' rhobarz(:,1)=rho(:,1) do l=2,nlay rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1)) enddo do l=1,nlay masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG enddo if (prt_level.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 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 !-------------------------------------------------------------------------------- ! if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out !===================================================================== ! Old version of thermcell_plume in thermcell_plume_6A.F90 ! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding ! to the 5B and 6A versions used for CMIP5 and CMIP6. ! The latest was previously named thermcellV1_plume. ! The new thermcell_plume is a clean version (removing obsolete ! options) of thermcell_plume_6A. ! The 3 versions are controled by ! flag_thermals_ed <= 9 thermcell_plume_6A ! <= 19 thermcell_plume_5B ! else thermcell_plume (default 20 for convergence with 6A) ! Fredho !===================================================================== if (iflag_thermals_ed<=9) then ! print*,'THERM NOUVELLE/NOUVELLE Arnaud' CALL thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & & ,lev_out,lunout1,igout) elseif (iflag_thermals_ed<=19) then ! print*,'THERM RIO et al 2010, version d Arnaud' CALL thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & & ,lev_out,lunout1,igout) else CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & & ,lev_out,lunout1,igout) endif if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_plum lmix ') if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume' if (prt_level.ge.10) then write(lunout1,*) 'Dans thermcell_main 2' write(lunout1,*) 'lmin ',lmin(igout) write(lunout1,*) 'lalim ',lalim(igout) write(lunout1,*) ' ig l alim_star entr_star detr_star f_star ' write(lunout1,'(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,lcong,lintercong,lmix,zw2, & & zlev,lmax,zmax,zmax0,zmix,wmax,zcong) ! Attention, w2 est transforme en sa racine carree dans cette routine ! Le probleme vient du fait que linter et lmix sont souvent egaux a 1. wmax_tmp=0. do l=1,nlay wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l)) enddo ! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmax ') if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height' !------------------------------------------------------------------------------- ! Fermeture,determination de f !------------------------------------------------------------------------------- ! ! CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & & lalim,lmin,zmax_sec,wmax_sec) call test_ltherm(ngrid,nlay,pplay,lmin,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry lalim ') if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry' if (prt_level.ge.10) then write(lunout1,*) 'Dans thermcell_main 1b' write(lunout1,*) 'lmin ',lmin(igout) write(lunout1,*) 'lalim ',lalim(igout) write(lunout1,*) ' ig l alim_star entr_star detr_star f_star ' write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) & & ,l=1,lalim(igout)+4) endif ! Choix de la fonction d'alimentation utilisee pour la fermeture. ! Apparemment sans importance alim_star_clos(:,:)=alim_star(:,:) alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:) ! !CR Appel de la fermeture seche if (iflag_thermals_closure.eq.1) then CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, & & zlev,lalim,alim_star_clos,zmax_sec,wmax_sec,f) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Appel avec les zmax et wmax tenant compte de la condensation ! Semble moins bien marcher else if (iflag_thermals_closure.eq.2) then CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, & & zlev,lalim,alim_star,zmax,wmax,f) endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure' if (tau_thermals>1.) then lambda=exp(-ptimestep/tau_thermals) f0=(1.-lambda)*f+lambda*f0 else f0=f endif ! Test valable seulement en 1D mais pas genant if (.not. (f0(1).ge.0.) ) then abort_message = '.not. (f0(1).ge.0.)' CALL abort_physic (modname,abort_message,1) endif !------------------------------------------------------------------------------- !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,lev_out,lunout1,igout) !IM 060508 & detr,zqla,zmax,lev_out,lunout,igout) if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux' call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lmax ') !------------------------------------------------------------------ ! On ne prend pas directement les profils issus des calculs precedents ! mais on s'autorise genereusement une relaxation vers ceci avec ! une constante de temps tau_thermals (typiquement 1800s). !------------------------------------------------------------------ if (tau_thermals>1.) then lambda=exp(-ptimestep/tau_thermals) fm0=(1.-lambda)*fm+lambda*fm0 entr0=(1.-lambda)*entr+lambda*entr0 detr0=(1.-lambda)*detr+lambda*detr0 else fm0=fm entr0=entr detr0=detr endif !------------------------------------------------------------------ ! 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 !c------------------------------------------------------------------ ! calcul du transport vertical !------------------------------------------------------------------ IF (iflag_thermals_down .GT. 0) THEN if (debut) print*,'WARNING !!! routine thermcell_down en cours de developpement' entrdn=fact_thermals_down*detr0 detrdn=fact_thermals_down*entr0 ! we want to transport potential temperature, total water and momentum CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zthl,zdthladj) CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,p_o,pdoadj) CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zu,pduadj) CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zv,pdvadj) ELSE !-------------------------------------------------------------- ! Temperature potentielle liquide effectivement mélangée par les thermiques do ll=1,nlay do ig=1,ngrid zthl(ig,ll)=ptemp(ig,ll)/zpspsk(ig,ll) enddo enddo call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & zthl,zdthladj,zta,lev_out) do ll=1,nlay do ig=1,ngrid z_o(ig,ll)=p_o(ig,ll) enddo enddo call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & z_o,pdoadj,z_oa,lev_out) #ifdef ISO ! C Risi: on utilise directement la meme routine do ixt=1,ntiso do ll=1,nlay DO ig=1,ngrid xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll) xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll) enddo enddo call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out) do ll=1,nlay DO ig=1,ngrid xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll) enddo enddo enddo #endif #ifdef ISO #ifdef ISOVERIF DO ll=1,nlay DO ig=1,ngrid if (iso_eau.gt.0) then call iso_verif_egalite(xtpo(iso_eau,ig,ll), & & p_o(ig,ll),'thermcell_main 594') call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), & & pdoadj(ig,ll),'thermcell_main 596') endif if (iso_HDO.gt.0) then call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) & & /p_o(ig,ll),'thermcell_main 610') endif enddo enddo !DO ll=1,nlay write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq' #endif #endif !------------------------------------------------------------------ ! calcul du transport vertical du moment horizontal !------------------------------------------------------------------ !IM 090508 if (dvdq == 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*dvdq,zmax & & ,fraca,zmax & & ,zu,zv,pduadj,pdvadj,zua,zva,lev_out) else ! calcul purement conservatif pour le transport de V call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse & & ,zu,pduadj,zua,lev_out) call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse & & ,zv,pdvadj,zva,lev_out) endif ENDIF ! print*,'13 OK convect8' do l=1,nlay do ig=1,ngrid pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) enddo enddo if (prt_level.ge.1) print*,'14 OK convect8' !------------------------------------------------------------------ ! Calculs de diagnostiques pour les sorties !------------------------------------------------------------------ !calcul de fraca pour les sorties if (sorties) then if (prt_level.ge.1) 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 ! WARNING !!! verifier que c'est bien ztemp_env qu'on veut là CHI=ztemp_env(ig,1)/(1669.0-122.0*z_o(ig,1)/zqsat(ig,1)-ztemp_env(ig,1)) pcon(ig)=pplay(ig,1)*(z_o(ig,1)/zqsat(ig,1))**CHI enddo !IM do k=1,nlay do k=1,nlay-1 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 !IM ierr=0 do ig=1,ngrid if (pcon(ig).le.pplay(ig,nlay)) then zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100. ierr=1 endif enddo ! if (ierr==1) then ! abort_message = 'thermcellV0_main: les thermiques vont trop haut ' ! CALL abort_physic (modname,abort_message,1) ! endif if (prt_level.ge.1) 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 if (prt_level.ge.1) 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 if (prt_level.ge.1) print*,'14d OK convect8' if (prt_level.ge.10)write(lunout,*) & & 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10' do l=1,nlay do ig=1,ngrid zf=fraca(ig,l) zf2=zf/(1.-zf) ! thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2 if(zw2(ig,l).gt.1.e-10) then wth2(ig,l)=zf2*(zw2(ig,l))**2 else wth2(ig,l)=0. endif 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.-p_o(ig,l)*1000.)**2 !test: on calcul q2/p_o=ratqsc ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(p_o(ig,l)*1000.)) enddo enddo !calcul des flux: q, thetal et thetav do l=1,nlay do ig=1,ngrid wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-p_o(ig,l)*1000.) wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l)) wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l)) enddo enddo !calcul du ratqscdiff if (prt_level.ge.1) print*,'14e OK convect8' var=0. vardiff=0. ratqsdiff(:,:)=0. do l=1,nlay do ig=1,ngrid if (l<=lalim(ig)) then var=var+alim_star(ig,l)*zqta(ig,l)*1000. endif enddo enddo if (prt_level.ge.1) print*,'14f OK convect8' do l=1,nlay do ig=1,ngrid if (l<=lalim(ig)) then zf=fraca(ig,l) zf2=zf/(1.-zf) vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2 endif enddo enddo if (prt_level.ge.1) print*,'14g OK convect8' do l=1,nlay do ig=1,ngrid ratqsdiff(ig,l)=sqrt(vardiff)/(p_o(ig,l)*1000.) enddo enddo endif if (prt_level.ge.1) print*,'thermcell_main FIN OK' !print*,'thermcell_main fin' RETURN end subroutine thermcell_main !============================================================================= !///////////////////////////////////////////////////////////////////////////// !============================================================================= subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,p_o,ztva, & ! in & zqla,f_star,zw2,comment) ! in !============================================================================= USE lmdz_thermcell_ini, ONLY: prt_level IMPLICIT NONE integer i, k, ngrid,nlay real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,p_o,ztva,zqla real, intent(in), dimension(ngrid,nlay) :: f_star,zw2 integer, intent(in), dimension(ngrid) :: long real seuil character*21 comment seuil=0.25 if (prt_level.ge.1) THEN print*,'WARNING !!! TEST ',comment endif return ! test sur la hauteur des thermiques ... do i=1,ngrid !IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then if (prt_level.ge.10) 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,nlay write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*p_o(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k) enddo endif enddo return end subroutine test_ltherm ! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP ! On transporte pbl_tke pour donner therm_tke ! Copie conforme de la subroutine DTKE dans physiq.F ecrite par Frederic Hourdin !======================================================================= !/////////////////////////////////////////////////////////////////////// !======================================================================= subroutine thermcell_tke_transport( & & ngrid,nlay,ptimestep,fm0,entr0,rg,pplev, & ! in & therm_tke_max) ! out USE lmdz_thermcell_ini, ONLY: prt_level implicit none !======================================================================= ! ! Calcul du transport verticale dans la couche limite en presence ! de "thermiques" explicitement representes ! calcul du dq/dt une fois qu'on connait les ascendances ! !======================================================================= integer ngrid,nlay real, intent(in) :: ptimestep real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev real, intent(in), dimension(ngrid,nlay) :: entr0 real, intent(in) :: rg real, intent(out), dimension(ngrid,nlay) :: therm_tke_max real detr0(ngrid,nlay) real masse0(ngrid,nlay) real masse(ngrid,nlay),fm(ngrid,nlay+1) real entr(ngrid,nlay) real q(ngrid,nlay) integer lev_out ! niveau pour les print real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1) integer ig,k lev_out=0 if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' ! calcul du detrainement do k=1,nlay detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k) masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG enddo ! Decalage vertical des entrainements et detrainements. masse(:,1)=0.5*masse0(:,1) entr(:,1)=0.5*entr0(:,1) detr(:,1)=0.5*detr0(:,1) fm(:,1)=0. do k=1,nlay-1 masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1)) entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1)) detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1)) fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k) enddo fm(:,nlay+1)=0. q(:,:)=therm_tke_max(:,:) !!! nrlmd le 16/09/2010 do ig=1,ngrid qa(ig,1)=q(ig,1) enddo !!! if (1==1) then do k=2,nlay do ig=1,ngrid if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. & & 1.e-5*masse(ig,k)) then qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & & /(fm(ig,k+1)+detr(ig,k)) else qa(ig,k)=q(ig,k) endif if (qa(ig,k).lt.0.) then ! print*,'qa<0!!!' endif if (q(ig,k).lt.0.) then ! print*,'q<0!!!' endif enddo enddo ! Calcul du flux subsident do k=2,nlay do ig=1,ngrid wqd(ig,k)=fm(ig,k)*q(ig,k) if (wqd(ig,k).lt.0.) then ! print*,'wqd<0!!!' endif enddo enddo do ig=1,ngrid wqd(ig,1)=0. wqd(ig,nlay+1)=0. enddo ! Calcul des tendances do k=1,nlay do ig=1,ngrid q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) & & -wqd(ig,k)+wqd(ig,k+1)) & & *ptimestep/masse(ig,k) enddo enddo endif therm_tke_max(:,:)=q(:,:) return !!! fin nrlmd le 10/04/2012 end subroutine thermcell_tke_transport END MODULE lmdz_thermcell_main