Index: /LMDZ6/trunk/libf/phylmd/alpale.h =================================================================== --- /LMDZ6/trunk/libf/phylmd/alpale.h (revision 4090) +++ /LMDZ6/trunk/libf/phylmd/alpale.h (revision 4090) @@ -0,0 +1,18 @@ +!===================================================================== +! Specifique de Ale/Alp : +!===================================================================== +! Mettre ca dans un alealp_th.sh en attendant mieux ... + +! dans alealp_th, thermcell_alp, physiq_mod, conf_phys + integer :: iflag_trig_bl,iflag_clos_bl + integer :: tau_trig_shallow,tau_trig_deep + real :: s_trig +! thermcell_alp et convection ... + integer :: iflag_coupl,iflag_clos,iflag_wake +! thermcell_alp + real :: alp_bl_k + + common/calpale1/iflag_trig_bl,iflag_clos_bl,tau_trig_shallow,tau_trig_deep + common/calpale2/s_trig,iflag_coupl,iflag_clos,iflag_wake,alp_bl_k + +!$OMP THREADPRIVATE(/calpale1/,/calpale2/) Index: DZ6/trunk/libf/phylmd/thermcell.h =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell.h (revision 4089) +++ (revision ) @@ -1,30 +1,0 @@ - integer :: iflag_thermals,nsplit_thermals - -!!! nrlmd le 10/04/2012 - integer :: iflag_trig_bl,iflag_clos_bl - integer :: tau_trig_shallow,tau_trig_deep - real :: s_trig -!!! fin nrlmd le 10/04/2012 - - real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30. - real :: alp_bl_k - real :: tau_thermals,fact_thermals_ed_dz - integer,parameter :: w2di_thermals=0 - integer :: isplit - - integer :: iflag_coupl,iflag_clos,iflag_wake - integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure - - common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure - common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz - common/ctherm4/iflag_coupl,iflag_clos,iflag_wake - common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux - -!!! nrlmd le 10/04/2012 - common/ctherm6/iflag_trig_bl,iflag_clos_bl - common/ctherm7/tau_trig_shallow,tau_trig_deep - common/ctherm8/s_trig -!!! fin nrlmd le 10/04/2012 - -!$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/) -!$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/) Index: DZ6/trunk/libf/phylmd/thermcellV0_main.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcellV0_main.F90 (revision 4089) +++ (revision ) @@ -1,2079 +1,0 @@ -! -! $Id$ -! - SUBROUTINE thermcellV0_main(itap,ngrid,nlay,ptimestep & - & ,pplay,pplev,pphi,debut & - & ,pu,pv,pt,po & - & ,pduadj,pdvadj,pdtadj,pdoadj & - & ,fm0,entr0,detr0,zqta,zqla,lmax & - & ,ratqscth,ratqsdiff,zqsatth & - & ,r_aspect,l_mix,tau_thermals & - & ,Ale_bl,Alp_bl,lalim_conv,wght_th & - & ,zmax0, f0,zw2,fraca) - - USE dimphy - USE print_control_mod, ONLY: prt_level,lunout - 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 "YOMCST.h" - include "YOETHF.h" - include "FCTTRE.h" - -! arguments: -! ---------- - -!IM 140508 - INTEGER itap - - INTEGER ngrid,nlay,w2di - real tau_thermals - 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 icount - data icount/0/ - save icount -!$OMP THREADPRIVATE(icount) - - integer,save :: igout=1 -!$OMP THREADPRIVATE(igout) - integer,save :: lunout1=6 -!$OMP THREADPRIVATE(lunout1) - integer,save :: lev_out=10 -!$OMP THREADPRIVATE(lev_out) - - INTEGER ig,k,l,ll - real zsortie1d(klon) - INTEGER lmax(klon),lmin(klon),lalim(klon) - INTEGER lmix(klon) - INTEGER lmix_bis(klon) - real linter(klon) - real zmix(klon) - real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1) -! real fraca(klon,klev) - - real zmax_sec(klon) -!on garde le zmax du pas de temps precedent - real zmax0(klon) -!FH/IM save zmax0 - - real lambda - - real zlev(klon,klev+1),zlay(klon,klev) - real deltaz(klon,klev) - REAL zh(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) -! FH probleme de dimensionnement avec l'allocation dynamique -! common/comtherm/thetath2,wth2 - - real ratqscth(klon,klev) - real var - real vardiff - real ratqsdiff(klon,klev) - - logical sorties - real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev) - real zpspsk(klon,klev) - - real wmax(klon) - real wmax_sec(klon) - real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev) - real fm(klon,klev+1),entr(klon,klev),detr(klon,klev) - - real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev) -!niveau de condensation - integer 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) -!FH/IM save f0 - real zlevinter(klon) - logical debut - real seuil - -! Declaration uniquement pour les sorties dans thermcell_out3d. -! Inutilise en 3D - real wthl(klon,klev) - real wthv(klon,klev) - real wq(klon,klev) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - -! - !nouvelles variables pour la convection - real Ale_bl(klon) - real Alp_bl(klon) - real alp_int(klon) - real ale_int(klon) - integer n_int(klon) - real fm_tot(klon) - real wght_th(klon,klev) - integer lalim_conv(klon) -!v1d logical therm -!v1d save therm - - character*2 str2 - character*10 str10 - - character (len=20) :: modname='thermcellV0_main' - character (len=80) :: abort_message - - EXTERNAL SCOPY -! - -!----------------------------------------------------------------------- -! initialisation: -! --------------- -! - - seuil=0.25 - - if (debut) then - fm0=0. - entr0=0. - detr0=0. - endif - fm=0. ; entr=0. ; detr=0. - - icount=icount+1 - -!IM 090508 beg -!print*,'=====================================================================' -!print*,'=====================================================================' -!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount -!print*,'=====================================================================' -!print*,'=====================================================================' -!IM 090508 end - - if (prt_level.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 -! - if (prt_level.ge.10)write(lunout,*) & - & 'WARNING thermcell_main f0=max(f0,1.e-2)' - do ig=1,klon - f0(ig)=max(f0(ig),1.e-2) - 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 (prt_level.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 - -!IM - 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 - -!calcul de la masse - 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 à 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 thermcellV0_init(ngrid,nlay,ztv,zlay,zlev, & - & lalim,lmin,alim_star,alim_star_tot,lev_out) - -call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ') -call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ') - - - if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init' - if (prt_level.ge.10) then - write(lunout1,*) 'Dans thermcell_main 1' - write(lunout1,*) 'lmin ',lmin(igout) - write(lunout1,*) 'lalim ',lalim(igout) - write(lunout1,*) ' ig l alim_star thetav' - write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) & - & ,ztv(igout,l),l=1,lalim(igout)+4) - endif - -!v1d do ig=1,klon -!v1d if (alim_star(ig,1).gt.1.e-10) then -!v1d therm=.true. -!v1d endif -!v1d enddo -!----------------------------------------------------------------------------- -! 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 thermcellV0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & - & lalim,lmin,zmax_sec,wmax_sec,lev_out) - -call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') -call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,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 - - - -!--------------------------------------------------------------------------------- -!calcul du melange et des variables dans le thermique -!-------------------------------------------------------------------------------- -! - if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out -!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & - CALL thermcellV0_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & - & zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, & - & lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, & - & ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter & - & ,lev_out,lunout1,igout) - if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out - - call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') - call testV0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,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,lmix,zw2, & - & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out) - - - call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') - call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') - call testV0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') - call testV0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ') - - if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height' - -!------------------------------------------------------------------------------- -! Fermeture,determination de f -!------------------------------------------------------------------------------- -! -!avant closure: on redéfinit lalim, alim_star_tot et alim_star -! do ig=1,klon -! do l=2,lalim(ig) -! alim_star(ig,l)=entr_star(ig,l) -! entr_star(ig,l)=0. -! enddo -! enddo - - CALL thermcellV0_closure(ngrid,nlay,r_aspect,ptimestep,rho, & - & zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out) - - 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 = 'Dans thermcell_main f0(1).lt.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 testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') - call testV0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,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 - -!c------------------------------------------------------------------ -! calcul du transport vertical -!------------------------------------------------------------------ - - call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse, & - & zthl,zdthladj,zta,lev_out) - call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse, & - & po,pdoadj,zoa,lev_out) - -!------------------------------------------------------------------ -! Calcul de la fraction de l'ascendance -!------------------------------------------------------------------ - do ig=1,klon - fraca(ig,1)=0. - fraca(ig,nlay+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 - -!------------------------------------------------------------------ -! calcul du transport vertical du moment horizontal -!------------------------------------------------------------------ - -!IM 090508 - if (1.eq.1) then -!IM 070508 vers. _dq -! 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) -!IM 050508 & ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out) - else - -! calcul purement conservatif pour le transport de V - call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse & - & ,zu,pduadj,zua,lev_out) - call thermcell_dq(ngrid,nlay,1,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 - - 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 - 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 -!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 - 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. - abort_message = 'thermcellV0_main: les thermiques vont trop haut ' - CALL abort_physic (modname,abort_message,1) - endif - enddo - 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*(zha(ig,l)-zh(ig,l)/zpspsk(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 -! 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 - - if (prt_level.ge.10) then - print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2 - ig=igout - do l=1,nlay - print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l) - enddo - do l=1,nlay - print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l) - enddo - endif - - do ig=1,ngrid - alp_int(ig)=0. - ale_int(ig)=0. - n_int(ig)=0 - enddo -! - do l=1,nlay - do ig=1,ngrid - if(l.LE.lmax(ig)) THEN - alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l) - ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2 - n_int(ig)=n_int(ig)+1 - endif - enddo - enddo -! print*,'avant calcul ale et alp' -!calcul de ALE et ALP pour la convection - do ig=1,ngrid -! Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig)) -! Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig)) -! Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig)) -! & *0.1 -!valeur integree de alp_bl * 0.5: - if (n_int(ig).gt.0) then - Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig) -! if (Alp_bl(ig).lt.0.) then -! Alp_bl(ig)=0. - endif -! endif -! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)), -! s wth3(ig,nivcon(ig)),Alp_bl(ig) -! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig) -! Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2 -! if (nivcon(ig).eq.1) then -! Ale_bl(ig)=0. -! else -!valeur max de ale_bl: - Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2 -! & /2. -! & *0.1 -! Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2 -! if (n_int(ig).gt.0) then -! Ale_bl(ig)=ale_int(ig)/n_int(ig) -! Ale_bl(ig)=4. -! endif -! endif -! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig)) -! Ale_bl(ig)=wth2(ig,nivcon(ig)) -! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig) - enddo -!test:calcul de la ponderation des couches pour KE -!initialisations -! print*,'ponderation' - do ig=1,ngrid - fm_tot(ig)=0. - enddo - do ig=1,ngrid - do k=1,klev - wght_th(ig,k)=1. - enddo - enddo - do ig=1,ngrid -! lalim_conv(ig)=lmix_bis(ig) -!la hauteur de la couche alim_conv = hauteur couche alim_therm - lalim_conv(ig)=lalim(ig) -! zentr(ig)=zlev(ig,lalim(ig)) - enddo - do ig=1,ngrid - do k=1,lalim_conv(ig) - fm_tot(ig)=fm_tot(ig)+fm(ig,k) - enddo - enddo - do ig=1,ngrid - do k=1,lalim_conv(ig) - if (fm_tot(ig).gt.1.e-10) then -! wght_th(ig,k)=fm(ig,k)/fm_tot(ig) - endif -!on pondere chaque couche par a* - if (alim_star(ig,k).gt.1.e-10) then - wght_th(ig,k)=alim_star(ig,k) - else - wght_th(ig,k)=1. - endif - enddo - enddo -! print*,'apres wght_th' -!test pour prolonger la convection - do ig=1,ngrid -!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then - if ((alim_star(ig,1).lt.1.e-10)) then - lalim_conv(ig)=1 - wght_th(ig,1)=1. -! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1) - endif - enddo - -!calcul du ratqscdiff - if (prt_level.ge.1) 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 - if (prt_level.ge.1) 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 - if (prt_level.ge.1) 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' - - if (prt_level.ge.1) print*,'thermcell_main sorties 3D' - endif - - if (prt_level.ge.1) print*,'thermcell_main FIN OK' - -! if(icount.eq.501) stop'au pas 301 dans thermcell_main' - return - end - -!----------------------------------------------------------------------------- - - subroutine testV0_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment) - USE print_control_mod, ONLY: prt_level - IMPLICIT NONE - - integer i, k, 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 - - if (prt_level.ge.1) THEN - print*,'WARNING !!! TEST ',comment - endif - return - -! test sur la hauteur des thermiques ... - do i=1,klon -!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,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 - endif - enddo - - - return - end - -!============================================================================== - SUBROUTINE thermcellV0_closure(ngrid,nlay,r_aspect,ptimestep,rho, & - & zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out) - -!------------------------------------------------------------------------- -!thermcell_closure: fermeture, determination de f -!------------------------------------------------------------------------- - USE print_control_mod, ONLY: prt_level,lunout - IMPLICIT NONE - - include "thermcell.h" - INTEGER ngrid,nlay - INTEGER ig,k - REAL r_aspect,ptimestep - integer lev_out ! niveau pour les print - - INTEGER lalim(ngrid) - REAL alim_star(ngrid,nlay) - REAL alim_star_tot(ngrid) - REAL rho(ngrid,nlay) - REAL zlev(ngrid,nlay) - REAL zmax(ngrid),zmax_sec(ngrid) - REAL wmax(ngrid),wmax_sec(ngrid) - real zdenom - - REAL alim_star2(ngrid) - - REAL f(ngrid) - - character (len=20) :: modname='thermcellV0_main' - character (len=80) :: abort_message - - do ig=1,ngrid - alim_star2(ig)=0. - enddo - do ig=1,ngrid - if (alim_star(ig,1).LT.1.e-10) then - f(ig)=0. - else - do k=1,lalim(ig) - alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2 & - & /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))) - enddo - zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig) - if (zdenom<1.e-14) then - print*,'ig=',ig - print*,'alim_star2',alim_star2(ig) - print*,'zmax',zmax(ig) - print*,'r_aspect',r_aspect - print*,'zdenom',zdenom - print*,'alim_star',alim_star(ig,:) - print*,'zmax_sec',zmax_sec(ig) - print*,'wmax_sec',wmax_sec(ig) - abort_message = 'zdenom<1.e-14' - CALL abort_physic (modname,abort_message,1) - endif - if ((zmax_sec(ig).gt.1.e-10).and.(iflag_thermals_ed.eq.0)) then - f(ig)=wmax_sec(ig)*alim_star_tot(ig)/(max(500.,zmax_sec(ig))*r_aspect & - & *alim_star2(ig)) -! f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ & -! & zmax_sec(ig))*wmax_sec(ig)) - if(prt_level.GE.10) write(lunout,*)'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig) - else - f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom -! f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ & -! & zmax(ig))*wmax(ig)) - if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig) - endif - endif -! f0(ig)=f(ig) - enddo - if (prt_level.ge.1) print*,'apres fermeture' - -! - return - end -!============================================================================== - SUBROUTINE thermcellV0_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, & - & zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, & - & lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, & - & ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter & - & ,lev_out,lunout1,igout) - -!-------------------------------------------------------------------------- -!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance -!-------------------------------------------------------------------------- - - USE print_control_mod, ONLY: prt_level - IMPLICIT NONE - - include "YOMCST.h" - include "YOETHF.h" - include "FCTTRE.h" - include "thermcell.h" - - INTEGER itap - INTEGER lunout1,igout - INTEGER ngrid,klev - REAL ptimestep - REAL ztv(ngrid,klev) - REAL zthl(ngrid,klev) - REAL po(ngrid,klev) - REAL zl(ngrid,klev) - REAL rhobarz(ngrid,klev) - REAL zlev(ngrid,klev+1) - REAL pplev(ngrid,klev+1) - REAL pphi(ngrid,klev) - REAL zpspsk(ngrid,klev) - REAL alim_star(ngrid,klev) - REAL zmax_sec(ngrid) - REAL f0(ngrid) - REAL l_mix - REAL r_aspect - INTEGER lalim(ngrid) - integer lev_out ! niveau pour les print - real zcon2(ngrid) - - real alim_star_tot(ngrid) - - REAL ztva(ngrid,klev) - REAL ztla(ngrid,klev) - REAL zqla(ngrid,klev) - REAL zqla0(ngrid,klev) - REAL zqta(ngrid,klev) - REAL zha(ngrid,klev) - - REAL detr_star(ngrid,klev) - REAL coefc - REAL detr_stara(ngrid,klev) - REAL detr_starb(ngrid,klev) - REAL detr_starc(ngrid,klev) - REAL detr_star0(ngrid,klev) - REAL detr_star1(ngrid,klev) - REAL detr_star2(ngrid,klev) - - REAL entr_star(ngrid,klev) - REAL entr_star1(ngrid,klev) - REAL entr_star2(ngrid,klev) - REAL detr(ngrid,klev) - REAL entr(ngrid,klev) - - REAL zw2(ngrid,klev+1) - REAL w_est(ngrid,klev+1) - REAL f_star(ngrid,klev+1) - REAL wa_moy(ngrid,klev+1) - - REAL ztva_est(ngrid,klev) - REAL zqla_est(ngrid,klev) - REAL zqsatth(ngrid,klev) - REAL zta_est(ngrid,klev) - - REAL linter(ngrid) - INTEGER lmix(ngrid) - INTEGER lmix_bis(ngrid) - REAL wmaxa(ngrid) - - INTEGER ig,l,k - - real zcor,zdelta,zcvm5,qlbef - real Tbef,qsatbef - real dqsat_dT,DT,num,denom - REAL REPS,RLvCp,DDT0 - PARAMETER (DDT0=.01) - logical Zsat - REAL fact_gamma,fact_epsilon - REAL c2(ngrid,klev) - - Zsat=.false. -! Initialisation - RLvCp = RLVTT/RCPD - - if (iflag_thermals_ed==0) then - fact_gamma=1. - fact_epsilon=1. - else if (iflag_thermals_ed==1) then - fact_gamma=1. - fact_epsilon=1. - else if (iflag_thermals_ed==2) then - fact_gamma=1. - fact_epsilon=2. - endif - - do l=1,klev - do ig=1,ngrid - zqla_est(ig,l)=0. - ztva_est(ig,l)=ztva(ig,l) - zqsatth(ig,l)=0. - enddo - enddo - -!CR: attention test couche alim -! do l=2,klev -! do ig=1,ngrid -! alim_star(ig,l)=0. -! enddo -! enddo -!AM:initialisations du thermique - do k=1,klev - do ig=1,ngrid - ztva(ig,k)=ztv(ig,k) - ztla(ig,k)=zthl(ig,k) - zqla(ig,k)=0. - zqta(ig,k)=po(ig,k) -! - ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k) - ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k) - zha(ig,k) = ztva(ig,k) -! - enddo - enddo - do k=1,klev - do ig=1,ngrid - detr_star(ig,k)=0. - entr_star(ig,k)=0. - - detr_stara(ig,k)=0. - detr_starb(ig,k)=0. - detr_starc(ig,k)=0. - detr_star0(ig,k)=0. - zqla0(ig,k)=0. - detr_star1(ig,k)=0. - detr_star2(ig,k)=0. - entr_star1(ig,k)=0. - entr_star2(ig,k)=0. - - detr(ig,k)=0. - entr(ig,k)=0. - enddo - enddo - if (prt_level.ge.1) print*,'7 OK convect8' - do k=1,klev+1 - do ig=1,ngrid - zw2(ig,k)=0. - w_est(ig,k)=0. - f_star(ig,k)=0. - wa_moy(ig,k)=0. - enddo - enddo - - if (prt_level.ge.1) print*,'8 OK convect8' - do ig=1,ngrid - linter(ig)=1. - lmix(ig)=1 - lmix_bis(ig)=2 - wmaxa(ig)=0. - enddo - -!----------------------------------------------------------------------------------- -!boucle de calcul de la vitesse verticale dans le thermique -!----------------------------------------------------------------------------------- - do l=1,klev-1 - do ig=1,ngrid - - - -! Calcul dans la premiere couche active du thermique (ce qu'on teste -! en disant que la couche est instable et que w2 en bas de la couche -! est nulle. - - if (ztv(ig,l).gt.ztv(ig,l+1) & - & .and.alim_star(ig,l).gt.1.e-10 & - & .and.zw2(ig,l).lt.1e-10) then - - -! Le panache va prendre au debut les caracteristiques de l'air contenu -! dans cette couche. - ztla(ig,l)=zthl(ig,l) - zqta(ig,l)=po(ig,l) - zqla(ig,l)=zl(ig,l) - f_star(ig,l+1)=alim_star(ig,l) - - zw2(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)) - w_est(ig,l+1)=zw2(ig,l+1) -! - - - else if ((zw2(ig,l).ge.1e-10).and. & - & (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then -!estimation du detrainement a partir de la geometrie du pas precedent -!tests sur la definition du detr -!calcul de detr_star et entr_star - - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! FH le test miraculeux de Catherine ? Le bout du tunel ? -! w_est(ig,3)=zw2(ig,2)* & -! & ((f_star(ig,2))**2) & -! & /(f_star(ig,2)+alim_star(ig,2))**2+ & -! & 2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) & -! & *(zlev(ig,3)-zlev(ig,2)) -! if (l.gt.2) then -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - -! Premier calcul de la vitesse verticale a partir de la temperature -! potentielle virtuelle - -! FH CESTQUOI CA ???? -#define int1d2 -!#undef int1d2 -#ifdef int1d2 - if (l.ge.2) then -#else - if (l.gt.2) then -#endif - - if (1.eq.1) then - w_est(ig,3)=zw2(ig,2)* & - & ((f_star(ig,2))**2) & - & /(f_star(ig,2)+alim_star(ig,2))**2+ & - & 2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) & -! & *1./3. & - & *(zlev(ig,3)-zlev(ig,2)) - endif - - -!--------------------------------------------------------------------------- -!calcul de l entrainement et du detrainement lateral -!--------------------------------------------------------------------------- -! -!test:estimation de ztva_new_est sans entrainement - - Tbef=ztla(ig,l-1)*zpspsk(ig,l) - zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) - qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l) - qsatbef=MIN(0.5,qsatbef) - zcor=1./(1.-retv*qsatbef) - qsatbef=qsatbef*zcor - Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10) - if (Zsat) then - qlbef=max(0.,zqta(ig,l-1)-qsatbef) - DT = 0.5*RLvCp*qlbef - do while (abs(DT).gt.DDT0) - Tbef=Tbef+DT - zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) - qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l) - qsatbef=MIN(0.5,qsatbef) - zcor=1./(1.-retv*qsatbef) - qsatbef=qsatbef*zcor - qlbef=zqta(ig,l-1)-qsatbef - - zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) - zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta - zcor=1./(1.-retv*qsatbef) - dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor) - num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef - denom=1.+RLvCp*dqsat_dT - DT=num/denom - enddo - zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef) - endif - ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) - ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) - zta_est(ig,l)=ztva_est(ig,l) - ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & - & -zqla_est(ig,l))-zqla_est(ig,l)) - - w_est(ig,l+1)=zw2(ig,l)* & - & ((f_star(ig,l))**2) & - & /(f_star(ig,l)+alim_star(ig,l))**2+ & - & 2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) & -! & *1./3. & - & *(zlev(ig,l+1)-zlev(ig,l)) - if (w_est(ig,l+1).lt.0.) then - w_est(ig,l+1)=zw2(ig,l) - endif -! -!calcul du detrainement -!======================= - -!CR:on vire les modifs - if (iflag_thermals_ed==0) then - -! Modifications du calcul du detrainement. -! Dans la version de la these de Catherine, on passe brusquement -! de la version seche a la version nuageuse pour le detrainement -! ce qui peut occasioner des oscillations. -! dans la nouvelle version, on commence par calculer un detrainement sec. -! Puis un autre en cas de nuages. -! Puis on combine les deux lineairement en fonction de la quantite d'eau. - -#define int1d3 -!#undef int1d3 -#define RIO_TH -#ifdef RIO_TH -!1. Cas non nuageux -! 1.1 on est sous le zmax_sec et w croit - if ((w_est(ig,l+1).gt.w_est(ig,l)).and. & - & (zlev(ig,l+1).lt.zmax_sec(ig)).and. & -#ifdef int1d3 - & (zqla_est(ig,l).lt.1.e-10)) then -#else - & (zqla(ig,l-1).lt.1.e-10)) then -#endif - detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) & - & *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1)) & - & -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l))) & - & /(r_aspect*zmax_sec(ig))) - detr_stara(ig,l)=detr_star(ig,l) - - if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l - -! 1.2 on est sous le zmax_sec et w decroit - else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and. & -#ifdef int1d3 - & (zqla_est(ig,l).lt.1.e-10)) then -#else - & (zqla(ig,l-1).lt.1.e-10)) then -#endif - detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) & - & /(rhobarz(ig,lmix(ig))*wmaxa(ig))* & - & (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) & - & *((zmax_sec(ig)-zlev(ig,l+1))/ & - & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. & - & -rhobarz(ig,l)*sqrt(w_est(ig,l)) & - & *((zmax_sec(ig)-zlev(ig,l))/ & - & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.) - detr_starb(ig,l)=detr_star(ig,l) - - if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l - - else - -! 1.3 dans les autres cas - detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l) & - & *(zlev(ig,l+1)-zlev(ig,l)) - detr_starc(ig,l)=detr_star(ig,l) - - if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l - - endif - -#else - -! 1.1 on est sous le zmax_sec et w croit - if ((w_est(ig,l+1).gt.w_est(ig,l)).and. & - & (zlev(ig,l+1).lt.zmax_sec(ig)) ) then - detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) & - & *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1)) & - & -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l))) & - & /(r_aspect*zmax_sec(ig))) - - if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l - -! 1.2 on est sous le zmax_sec et w decroit - else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then - detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) & - & /(rhobarz(ig,lmix(ig))*wmaxa(ig))* & - & (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) & - & *((zmax_sec(ig)-zlev(ig,l+1))/ & - & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. & - & -rhobarz(ig,l)*sqrt(w_est(ig,l)) & - & *((zmax_sec(ig)-zlev(ig,l))/ & - & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.) - if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l - - else - detr_star=0. - endif - -! 1.3 dans les autres cas - detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l) & - & *(zlev(ig,l+1)-zlev(ig,l)) - - coefc=min(zqla(ig,l-1)/1.e-3,1.) - if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1. - coefc=1. -! il semble qu'il soit important de baser le calcul sur -! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l) - detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc) - - if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l - -#endif - - - if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l -!IM 730508 beg -! if(itap.GE.7200) THEN -! print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l) -! endif -!IM 730508 end - - zqla0(ig,l)=zqla_est(ig,l) - detr_star0(ig,l)=detr_star(ig,l) -!IM 060508 beg -! if(detr_star(ig,l).GT.1.) THEN -! print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) & -! & ,detr_starc(ig,l),coefc -! endif -!IM 060508 end -!IM 160508 beg -!IM 160508 IF (f0(ig).NE.0.) THEN - detr_star(ig,l)=detr_star(ig,l)/f0(ig) -!IM 160508 ELSE IF(detr_star(ig,l).EQ.0.) THEN -!IM 160508 print*,'WARNING1 : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap -!IM 160508 ELSE -!IM 160508 print*,'WARNING2 : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l) -!IM 160508 ENDIF -!IM 160508 end -!IM 060508 beg -! if(detr_star(ig,l).GT.1.) THEN -! print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), & -! & REAL(1)/f0(ig) -! endif -!IM 060508 end - if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l -! -!calcul de entr_star - -! #undef test2 -! #ifdef test2 -! La version test2 destabilise beaucoup le modele. -! Il semble donc que ca aide d'avoir un entrainement important sous -! le nuage. -! if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then -! entr_star(ig,l)=0.4*detr_star(ig,l) -! else -! entr_star(ig,l)=0. -! endif -! #else -! -! Deplacement du calcul de entr_star pour eviter d'avoir aussi -! entr_star > fstar. -! Redeplacer suite a la transformation du cas detr>f -! FH - - if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l -#define int1d -!FH 070508 #define int1d4 -!#undef int1d4 -! L'option int1d4 correspond au choix dans le cas ou le detrainement -! devient trop grand. - -#ifdef int1d - -#ifdef int1d4 -#else - detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l)) -!FH 070508 plus - detr_star(ig,l)=min(detr_star(ig,l),1.) -#endif - - entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.) - - if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l -#ifdef int1d4 -! Si le detrainement excede le flux en bas + l'entrainement, le thermique -! doit disparaitre. - if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then - detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l) - f_star(ig,l+1)=0. - linter(ig)=l+1 - zw2(ig,l+1)=-1.e-10 - endif -#endif - - -#else - - if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l - if(l.gt.lalim(ig)) then - entr_star(ig,l)=0.4*detr_star(ig,l) - else - -! FH : -! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1 -! en haut de la couche d'alimentation. -! A remettre en questoin a la premiere occasion mais ca peut aider a -! ecrire un code robuste. -! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais -! d* non nul) on a une discontinuité de e* ou d* en haut de la couche -! d'alimentation, ce qui n'est pas forcement heureux. - - if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l -#undef pre_int1c -#ifdef pre_int1c - entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.) - detr_star(ig,l)=entr_star(ig,l) -#else - entr_star(ig,l)=0. -#endif - - endif - -#endif - - if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l - entr_star1(ig,l)=entr_star(ig,l) - detr_star1(ig,l)=detr_star(ig,l) -! - -#ifdef int1d -#else - if (detr_star(ig,l).gt.f_star(ig,l)) then - -! Ce test est là entre autres parce qu'on passe par des valeurs -! delirantes de detr_star. -! ca vaut sans doute le coup de verifier pourquoi. - - detr_star(ig,l)=f_star(ig,l) -#ifdef pre_int1c - if (l.gt.lalim(ig)+1) then - entr_star(ig,l)=0. - alim_star(ig,l)=0. -! FH ajout pour forcer a stoper le thermique juste sous le sommet -! de la couche (voir calcul de finter) - zw2(ig,l+1)=-1.e-10 - linter(ig)=l+1 - else - entr_star(ig,l)=0.4*detr_star(ig,l) - endif -#else - entr_star(ig,l)=0.4*detr_star(ig,l) -#endif - endif -#endif - - else !l > 2 - detr_star(ig,l)=0. - entr_star(ig,l)=0. - endif - - entr_star2(ig,l)=entr_star(ig,l) - detr_star2(ig,l)=detr_star(ig,l) - if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l - - endif ! iflag_thermals_ed==0 - -!CR:nvlle def de entr_star et detr_star - if (iflag_thermals_ed>=1) then -! if (l.lt.lalim(ig)) then -! if (l.lt.2) then -! entr_star(ig,l)=0. -! detr_star(ig,l)=0. -! else -! if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then -! entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) -! else -! entr_star(ig,l)= & -! & f_star(ig,l)/(2.*w_est(ig,l+1)) & -! & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) & -! & *(zlev(ig,l+1)-zlev(ig,l)) - - - entr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)), & - & f_star(ig,l)/(2.*w_est(ig,l+1)) & - & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) & - & *(zlev(ig,l+1)-zlev(ig,l))) & - & +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) - - if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then - alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l) - lalim(ig)=lmix_bis(ig) - if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l) - endif - - if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then -! c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l)) - c2(ig,l)=0.001 - detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)), & - & c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) & - & -f_star(ig,l)/(2.*w_est(ig,l+1)) & - & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) & - & *(zlev(ig,l+1)-zlev(ig,l))) & - & +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) - - else -! c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l)) - c2(ig,l)=0.003 - - detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)), & - & c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) & - & -f_star(ig,l)/(2.*w_est(ig,l+1)) & - & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) & - & *(zlev(ig,l+1)-zlev(ig,l))) & - & +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) - endif - - -! detr_star(ig,l)=detr_star(ig,l)*3. -! if (l.lt.lalim(ig)) then -! entr_star(ig,l)=0. -! endif -! if (l.lt.2) then -! entr_star(ig,l)=0. -! detr_star(ig,l)=0. -! endif - - -! endif -! else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then -! entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1)) & -! & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) & -! & *(zlev(ig,l+1)-zlev(ig,l)) -! detr_star(ig,l)=0.002*f_star(ig,l) & -! & *(zlev(ig,l+1)-zlev(ig,l)) -! else -! entr_star(ig,l)=0.001*f_star(ig,l) & -! & *(zlev(ig,l+1)-zlev(ig,l)) -! detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1)) & -! & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) & -! & *(zlev(ig,l+1)-zlev(ig,l)) & -! & +0.002*f_star(ig,l) & -! & *(zlev(ig,l+1)-zlev(ig,l)) -! endif - - endif ! iflag_thermals_ed==1 - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! FH inutile si on conserve comme on l'a fait plus haut entr=detr -! dans la couche d'alimentation -!pas d entrainement dans la couche alim -! if ((l.le.lalim(ig))) then -! entr_star(ig,l)=0. -! endif -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! -!prise en compte du detrainement et de l entrainement dans le calcul du flux - - f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & - & -detr_star(ig,l) - -!test sur le signe de f_star - if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l - if (f_star(ig,l+1).gt.1.e-10) then -!---------------------------------------------------------------------------- -!calcul de la vitesse verticale en melangeant Tl et qt du thermique -!--------------------------------------------------------------------------- -! - Zsat=.false. - ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & - & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & - & /(f_star(ig,l+1)+detr_star(ig,l)) -! - zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & - & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & - & /(f_star(ig,l+1)+detr_star(ig,l)) -! - Tbef=ztla(ig,l)*zpspsk(ig,l) - zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) - qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l) - qsatbef=MIN(0.5,qsatbef) - zcor=1./(1.-retv*qsatbef) - qsatbef=qsatbef*zcor - Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10) - if (Zsat) then - qlbef=max(0.,zqta(ig,l)-qsatbef) - DT = 0.5*RLvCp*qlbef - do while (abs(DT).gt.DDT0) - Tbef=Tbef+DT - zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) - qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l) - qsatbef=MIN(0.5,qsatbef) - zcor=1./(1.-retv*qsatbef) - qsatbef=qsatbef*zcor - qlbef=zqta(ig,l)-qsatbef - - zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) - zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta - zcor=1./(1.-retv*qsatbef) - dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor) - num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef - denom=1.+RLvCp*dqsat_dT - DT=num/denom - enddo - zqla(ig,l) = max(0.,qlbef) - endif -! - if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l -! on ecrit de maniere conservative (sat ou non) -! T = Tl +Lv/Cp ql - ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) - ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) -!on rajoute le calcul de zha pour diagnostiques (temp potentielle) - zha(ig,l) = ztva(ig,l) - ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & - & -zqla(ig,l))-zqla(ig,l)) - -!on ecrit zqsat - zqsatth(ig,l)=qsatbef -!calcul de vitesse - zw2(ig,l+1)=zw2(ig,l)* & - & ((f_star(ig,l))**2) & -! Tests de Catherine -! & /(f_star(ig,l+1)+detr_star(ig,l))**2+ & - & /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ & - & 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) & - & *fact_gamma & - & *(zlev(ig,l+1)-zlev(ig,l)) -!prise en compte des forces de pression que qd flottabilité<0 -! zw2(ig,l+1)=zw2(ig,l)* & -! & 1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + & -! & (f_star(ig,l))**2 & -! & /(f_star(ig,l)+entr_star(ig,l))**2+ & -! & (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+ & -! & /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ & -! & /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ & -! & 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) & -! & *1./3. & -! & *(zlev(ig,l+1)-zlev(ig,l)) - -! write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), & -! & -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), & -! & 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) - - -! zw2(ig,l+1)=zw2(ig,l)* & -! & (2.-2.*entr_star(ig,l)/f_star(ig,l)) & -! & -zw2(ig,l-1)+ & -! & 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) & -! & *1./3. & -! & *(zlev(ig,l+1)-zlev(ig,l)) - - endif - endif - if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l -! -!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max - - 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).gt.0.).and. (zw2(ig,l+1).le.0.)) then - 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 - if (zqla(ig,l).lt.1.e-10) then - lmix_bis(ig)=l+1 - endif - lmix(ig)=l+1 - wmaxa(ig)=wa_moy(ig,l+1) - endif - enddo - enddo - -!on remplace a* par e* ds premiere couche -! if (iflag_thermals_ed.ge.1) then -! do ig=1,ngrid -! do l=2,klev -! if (l.lt.lalim(ig)) then -! alim_star(ig,l)=entr_star(ig,l) -! endif -! enddo -! enddo -! do ig=1,ngrid -! lalim(ig)=lmix_bis(ig) -! enddo -! endif - if (iflag_thermals_ed.ge.1) then - do ig=1,ngrid - do l=2,lalim(ig) - alim_star(ig,l)=entr_star(ig,l) - entr_star(ig,l)=0. - enddo - enddo - endif - if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l - -! print*,'thermcell_plume OK' - - return - end -!============================================================================== - SUBROUTINE thermcellV0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & - & lalim,lmin,zmax,wmax,lev_out) - -!-------------------------------------------------------------------------- -!thermcell_dry: calcul de zmax et wmax du thermique sec -!-------------------------------------------------------------------------- - USE print_control_mod, ONLY: prt_level - IMPLICIT NONE - include "YOMCST.h" - INTEGER l,ig - - INTEGER ngrid,nlay - REAL zlev(ngrid,nlay+1) - REAL pphi(ngrid,nlay) - REAl ztv(ngrid,nlay) - REAL alim_star(ngrid,nlay) - INTEGER lalim(ngrid) - integer lev_out ! niveau pour les print - - REAL zmax(ngrid) - REAL wmax(ngrid) - -!variables locales - REAL zw2(ngrid,nlay+1) - REAL f_star(ngrid,nlay+1) - REAL ztva(ngrid,nlay+1) - REAL wmaxa(ngrid) - REAL wa_moy(ngrid,nlay+1) - REAL linter(ngrid),zlevinter(ngrid) - INTEGER lmix(ngrid),lmax(ngrid),lmin(ngrid) - -!initialisations - do ig=1,ngrid - do l=1,nlay+1 - zw2(ig,l)=0. - wa_moy(ig,l)=0. - enddo - enddo - do ig=1,ngrid - do l=1,nlay - ztva(ig,l)=ztv(ig,l) - enddo - enddo - do ig=1,ngrid - wmax(ig)=0. - wmaxa(ig)=0. - enddo -!calcul de la vitesse a partir de la CAPE en melangeant thetav - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! A eliminer -! Ce if complique etait fait pour reperer la premiere couche instable -! Ici, c'est lmin. -! -! do l=1,nlay-2 -! do ig=1,ngrid -! if (ztv(ig,l).gt.ztv(ig,l+1) & -! & .and.alim_star(ig,l).gt.1.e-10 & -! & .and.zw2(ig,l).lt.1e-10) then -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - -! Calcul des F^*, integrale verticale de E^* - f_star(:,1)=0. - do l=1,nlay - f_star(:,l+1)=f_star(:,l)+alim_star(:,l) - enddo - -! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise - linter(:)=0. - -! couche la plus haute concernee par le thermique. - lmax(:)=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(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) -!------------------------------------------------------------------------ - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! A eliminer : dans cette version, si zw2 est > 0 on a un therique. -! et donc, au dessus, f_star(ig,l+1) est forcement suffisamment -! grand puisque on n'a pas de detrainement. -! f_star est une fonction croissante. -! c'est donc vraiment sur zw2 uniquement qu'il faut faire le test. -! else if ((zw2(ig,l).ge.1e-10).and. & -! & (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then -! f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - else if (zw2(ig,l).ge.1e-10) then - - ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+alim_star(ig,l) & - & *ztv(ig,l))/f_star(ig,l+1) - zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+ & - & 2.*RG*(ztva(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(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then -! print*,'On tombe sur le cas particulier de thermcell_dry' - zw2(ig,l+1)=0. - linter(ig)=l+1 - lmax(ig)=l - 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. - lmax(ig)=l - 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 - lmix(ig)=l+1 - wmaxa(ig)=wa_moy(ig,l+1) - endif - enddo - enddo - if (prt_level.ge.1) print*,'fin calcul zw2' -! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! A eliminer : -! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut -! Calcul de la couche correspondant a la hauteur 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 -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -! -! 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 - 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 - do ig=1,ngrid -! calcul de zlevinter - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! FH A eliminer -! Simplification -! zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))* & -! & linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) & -! & -zlev(ig,lmax(ig))) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - zlevinter(ig)=zlev(ig,lmax(ig)) + & - & (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig))) - zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig))) - enddo - -! Verification que lalim<=lmax - do ig=1,ngrid - if(lalim(ig)>lmax(ig)) then - if ( prt_level > 1 ) THEN - print*,'WARNING thermcell_dry ig=',ig,' lalim=',lalim(ig),' lmax(ig)=',lmax(ig) - endif - lmax(ig)=lalim(ig) - endif - enddo - - RETURN - END -!============================================================================== - SUBROUTINE thermcellV0_init(ngrid,nlay,ztv,zlay,zlev, & - & lalim,lmin,alim_star,alim_star_tot,lev_out) - -!---------------------------------------------------------------------- -!thermcell_init: calcul du profil d alimentation du thermique -!---------------------------------------------------------------------- - USE print_control_mod, ONLY: prt_level - IMPLICIT NONE - include "thermcell.h" - - INTEGER l,ig -!arguments d entree - INTEGER ngrid,nlay - REAL ztv(ngrid,nlay) - REAL zlay(ngrid,nlay) - REAL zlev(ngrid,nlay+1) -!arguments de sortie - INTEGER lalim(ngrid) - INTEGER lmin(ngrid) - REAL alim_star(ngrid,nlay) - REAL alim_star_tot(ngrid) - integer lev_out ! niveau pour les print - - REAL zzalim(ngrid) -!CR: ponderation entrainement des couches instables -!def des alim_star tels que alim=f*alim_star - - do l=1,nlay - do ig=1,ngrid - alim_star(ig,l)=0. - enddo - enddo -! determination de la longueur de la couche d entrainement - do ig=1,ngrid - lalim(ig)=1 - enddo - - if (iflag_thermals_ed.ge.1) then -!si la première couche est instable, on declenche un thermique - do ig=1,ngrid - if (ztv(ig,1).gt.ztv(ig,2)) then - lmin(ig)=1 - lalim(ig)=2 - alim_star(ig,1)=1. - alim_star_tot(ig)=alim_star(ig,1) - if(prt_level.GE.10) print*,'init',alim_star(ig,1),alim_star_tot(ig) - else - lmin(ig)=1 - lalim(ig)=1 - alim_star(ig,1)=0. - alim_star_tot(ig)=0. - endif - enddo - - else -!else iflag_thermals_ed=0 ancienne def de l alim - -!on ne considere que les premieres couches instables - do l=nlay-2,1,-1 - do ig=1,ngrid - if (ztv(ig,l).gt.ztv(ig,l+1).and. & - & ztv(ig,l+1).le.ztv(ig,l+2)) then - lalim(ig)=l+1 - endif - enddo - enddo - -! determination du lmin: couche d ou provient le thermique - - do ig=1,ngrid -! FH initialisation de lmin a nlay plutot que 1. -! lmin(ig)=nlay - lmin(ig)=1 - enddo - do l=nlay,2,-1 - do ig=1,ngrid - if (ztv(ig,l-1).gt.ztv(ig,l)) then - lmin(ig)=l-1 - endif - enddo - enddo -! - zzalim(:)=0. - do l=1,nlay-1 - do ig=1,ngrid - if (l1) then - zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig))) - else - zzalim(ig)=zlay(ig,1) - endif - enddo - - if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1)) - -! definition de l'entrainement des couches - if (1.eq.1) then - do l=1,nlay-1 - do ig=1,ngrid - if (ztv(ig,l).gt.ztv(ig,l+1).and. & - & l.ge.lmin(ig).and.l.lt.lalim(ig)) then -!def possibles pour alim_star: zdthetadz, dthetadz, zdtheta - alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & - & *sqrt(zlev(ig,l+1)) - endif - enddo - enddo - else - do l=1,nlay-1 - do ig=1,ngrid - if (ztv(ig,l).gt.ztv(ig,l+1).and. & - & l.ge.lmin(ig).and.l.lt.lalim(ig)) then - alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) & - & *(zlev(ig,l+1)-zlev(ig,l)) - endif - enddo - enddo - endif - -! pas de thermique si couche 1 stable - do ig=1,ngrid -!CRnouveau test - if (alim_star(ig,1).lt.1.e-10) then - do l=1,nlay - alim_star(ig,l)=0. - enddo - lmin(ig)=1 - endif - enddo -! calcul de l alimentation totale - do ig=1,ngrid - alim_star_tot(ig)=0. - enddo - do l=1,nlay - do ig=1,ngrid - alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) - enddo - enddo -! -! Calcul entrainement normalise - do l=1,nlay - do ig=1,ngrid - if (alim_star_tot(ig).gt.1.e-10) then - alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) - endif - enddo - enddo - -!on remet alim_star_tot a 1 - do ig=1,ngrid - alim_star_tot(ig)=1. - enddo - - endif -!endif iflag_thermals_ed - return - end -!============================================================================== Index: DZ6/trunk/libf/phylmd/thermcell_condens.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_condens.F90 (revision 4089) +++ (revision ) @@ -1,83 +1,0 @@ -subroutine thermcell_condens(klon,active,zpspsk,pplev,ztla,zqta,zqla) -implicit none - -#include "YOMCST.h" -#include "YOETHF.h" -#include "FCTTRE.h" - - -!==================================================================== -! DECLARATIONS -!==================================================================== - -! Arguments -INTEGER klon -REAL zpspsk(klon),pplev(klon) -REAL ztla(klon),zqta(klon),zqla(klon) -LOGICAL active(klon) - -! Variables locales -INTEGER ig,iter -REAL Tbef(klon),DT(klon) -REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT -logical Zsat -REAL RLvCp -REAL, SAVE :: DDT0=.01 - !$OMP THREADPRIVATE(DDT0) - -LOGICAL afaire(klon),tout_converge - -!==================================================================== -! INITIALISATIONS -!==================================================================== - -RLvCp = RLVTT/RCPD -tout_converge=.false. -afaire(:)=.false. -DT(:)=0. - - -!==================================================================== -! Routine a vectoriser en copiant active dans converge et en mettant -! la boucle sur les iterations a l'exterieur est en mettant -! converge= false des que la convergence est atteinte. -!==================================================================== - -do ig=1,klon - if (active(ig)) then - Tbef(ig)=ztla(ig)*zpspsk(ig) - zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig))) - qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig) - qsatbef=MIN(0.5,qsatbef) - zcor=1./(1.-retv*qsatbef) - qsatbef=qsatbef*zcor - qlbef=max(0.,zqta(ig)-qsatbef) - DT(ig) = 0.5*RLvCp*qlbef - endif -enddo - -do iter=1,10 - afaire(:)=abs(DT(:)).gt.DDT0 - do ig=1,klon - if (afaire(ig)) then - Tbef(ig)=Tbef(ig)+DT(ig) - zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig))) - qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig) - qsatbef=MIN(0.5,qsatbef) - zcor=1./(1.-retv*qsatbef) - qsatbef=qsatbef*zcor - qlbef=zqta(ig)-qsatbef - zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig))) - zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta - zcor=1./(1.-retv*qsatbef) - dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor) - num=-Tbef(ig)+ztla(ig)*zpspsk(ig)+RLvCp*qlbef - denom=1.+RLvCp*dqsat_dT - zqla(ig) = max(0.,zqta(ig)-qsatbef) - DT(ig)=num/denom - endif - enddo -enddo - -return -end Index: DZ6/trunk/libf/phylmd/thermcell_flux.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_flux.F90 (revision 4089) +++ (revision ) @@ -1,519 +1,0 @@ -! -! $Id$ -! - - - SUBROUTINE thermcell_flux(ngrid,klev,ptimestep,masse, & - & lalim,lmax,alim_star, & - & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, & - & detr,zqla,zmax,lev_out,lunout1,igout) - - -!--------------------------------------------------------------------------- -!thermcell_flux: deduction des flux -!--------------------------------------------------------------------------- - - USE print_control_mod, ONLY: prt_level,lunout - IMPLICIT NONE - - INTEGER ig,l - INTEGER ngrid,klev - - REAL alim_star(ngrid,klev),entr_star(ngrid,klev) - REAL detr_star(ngrid,klev) - REAL zw2(ngrid,klev+1) - REAL zlev(ngrid,klev+1) - REAL masse(ngrid,klev) - REAL ptimestep - REAL rhobarz(ngrid,klev) - REAL f(ngrid) - INTEGER lmax(ngrid) - INTEGER lalim(ngrid) - REAL zqla(ngrid,klev) - REAL zmax(ngrid) - - integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha - integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8 - - - REAL entr(ngrid,klev),detr(ngrid,klev) - REAL fm(ngrid,klev+1) - REAL zfm - - integer igout - integer lev_out - integer lunout1 - - REAL f_old,ddd0,eee0,ddd,eee,zzz - - REAL fomass_max,alphamax - save fomass_max,alphamax -!$OMP THREADPRIVATE(fomass_max,alphamax) - - character (len=20) :: modname='thermcell_flux' - character (len=80) :: abort_message - - fomass_max=0.5 - alphamax=0.7 - - ncorecfm1=0 - ncorecfm2=0 - ncorecfm3=0 - ncorecfm4=0 - ncorecfm5=0 - ncorecfm6=0 - ncorecfm7=0 - ncorecfm8=0 - ncorecalpha=0 - -!initialisation - fm(:,:)=0. - - if (prt_level.ge.10) then - write(lunout,*) 'Dans thermcell_flux 0' - write(lunout,*) 'flux base ',f(igout) - write(lunout,*) 'lmax ',lmax(igout) - write(lunout,*) 'lalim ',lalim(igout) - write(lunout,*) 'ig= ',igout - write(lunout,*) ' l E* A* D* ' - write(lunout,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) & - & ,l=1,lmax(igout)) - endif - - -!------------------------------------------------------------------------- -! Verification de la nullite des entrainement et detrainement au dessus -! de lmax(ig) -!------------------------------------------------------------------------- - if ( prt_level > 1 ) THEN - do l=1,klev - do ig=1,ngrid - if (l.le.lmax(ig)) then - if (entr_star(ig,l).gt.1.) then - print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig) - print*,'entr_star(ig,l)',entr_star(ig,l) - print*,'alim_star(ig,l)',alim_star(ig,l) - print*,'detr_star(ig,l)',detr_star(ig,l) - endif - else - if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then - print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig) - print*,'entr_star(ig,l)',entr_star(ig,l) - print*,'alim_star(ig,l)',alim_star(ig,l) - print*,'detr_star(ig,l)',detr_star(ig,l) - abort_message = '' - CALL abort_physic (modname,abort_message,1) - endif - endif - enddo - enddo - endif !( prt_level > 1 ) THEN -!------------------------------------------------------------------------- -! Multiplication par le flux de masse issu de la femreture -!------------------------------------------------------------------------- - - do l=1,klev - entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l)) - detr(:,l)=f(:)*detr_star(:,l) - enddo - - if (prt_level.ge.10) then - write(lunout,*) 'Dans thermcell_flux 1' - write(lunout,*) 'flux base ',f(igout) - write(lunout,*) 'lmax ',lmax(igout) - write(lunout,*) 'lalim ',lalim(igout) - write(lunout,*) 'ig= ',igout - write(lunout,*) ' l E D W2' - write(lunout,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) & - & ,zw2(igout,l+1),l=1,lmax(igout)) - endif - - fm(:,1)=0. - do l=1,klev - 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,klev - do ig=1,ngrid - if (detr(ig,l).gt.fm(ig,l)) then - ncorecfm8=ncorecfm8+1 -! igout=ig - endif - enddo - enddo - - if (prt_level.ge.10) & - & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & - & ptimestep,masse,entr,detr,fm,'2 ') - - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! 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,klev - - 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 l=1,klev - 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 -! enddo - - if (prt_level.ge.10) & - & write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & - & entr(igout,l),detr(igout,l),fm(igout,l+1) - -!------------------------------------------------------------------------- -!Test sur fraca croissant -!------------------------------------------------------------------------- - - - if (1.eq.1) then -! do l=1,klev - 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 -! enddo - - if (prt_level.ge.10) & - & write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & - & entr(igout,l),detr(igout,l),fm(igout,l+1) - else - if (l.eq.1) then - print*,'Test sur les fractions croissantes inhibe dans thermcell_flux2' - endif - endif - - -!------------------------------------------------------------------------- -!test sur flux de masse croissant -!------------------------------------------------------------------------- - -! do l=1,klev - 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 -! enddo - - if (prt_level.ge.10) & - & write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & - & entr(igout,l),detr(igout,l),fm(igout,l+1) - -!------------------------------------------------------------------------- -!detr ne peut pas etre superieur a fm -!------------------------------------------------------------------------- - - if(1.eq.1) then - -! do l=1,klev - do ig=1,ngrid - if (entr(ig,l)<0.) then - print*,'N1 ig,l,entr',ig,l,entr(ig,l) - abort_message = 'entr negatif' - CALL abort_physic (modname,abort_message,1) - endif - 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. - if (l.gt.lalim(ig)) then - lmax(ig)=l - fm(ig,l+1)=0. - entr(ig,l)=0. - else - ncorecfm7=ncorecfm7+1 - endif - endif - - if(l.gt.lmax(ig)) then - detr(ig,l)=0. - fm(ig,l+1)=0. - entr(ig,l)=0. - endif - - if (entr(ig,l).lt.0.) then - print*,'ig,l,lmax(ig)',ig,l,lmax(ig) - print*,'entr(ig,l)',entr(ig,l) - print*,'fm(ig,l)',fm(ig,l) - abort_message = 'probleme dans thermcell flux' - CALL abort_physic (modname,abort_message,1) - endif - enddo -! enddo - endif - - - if (prt_level.ge.10) & - & write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & - & entr(igout,l),detr(igout,l),fm(igout,l+1) - -!------------------------------------------------------------------------- -!fm ne peut pas etre negatif -!------------------------------------------------------------------------- - -! do l=1,klev - 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. -! print*,'fm2<0',l+1,lmax(ig) - ncorecfm2=ncorecfm2+1 - endif - if (detr(ig,l).lt.0.) then - print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig) - print*,'detr(ig,l)',detr(ig,l) - print*,'fm(ig,l)',fm(ig,l) - abort_message = 'probleme dans thermcell flux' - CALL abort_physic (modname,abort_message,1) - endif - enddo -! enddo - - if (prt_level.ge.10) & - & write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & - & entr(igout,l),detr(igout,l),fm(igout,l+1) - -!----------------------------------------------------------------------- -!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 l=1,klev - 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 -! zw2(ig,l+1)=0. -! zqla(ig,l+1)=0. - detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) -! lmax(ig)=l+1 -! zmax(ig)=zlev(ig,lmax(ig)) -! print*,'alpha>1',l+1,lmax(ig) - ncorecalpha=ncorecalpha+1 - endif - endif - enddo -! enddo -! - - - if (prt_level.ge.10) & - & write(lunout,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, & - & entr(igout,l),detr(igout,l),fm(igout,l+1) - -! Fin de la grande boucle sur les niveaux verticaux - enddo - - if (prt_level.ge.10) & - & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & - & ptimestep,masse,entr,detr,fm,'8 ') - - -!----------------------------------------------------------------------- -! 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 -!----------------------------------------------------------------------- - - if (1.eq.1) then - do l=1,klev-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 - if(l.ge.lmax(ig).and.0.eq.1) then - print*,'ig,l',ig,l - print*,'eee0',eee0 - print*,'ddd0',ddd0 - print*,'eee',eee - print*,'ddd',ddd - print*,'entr',entr(ig,l) - print*,'detr',detr(ig,l) - print*,'masse',masse(ig,l) - print*,'fomass_max',fomass_max - print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep - print*,'ptimestep',ptimestep - print*,'lmax(ig)',lmax(ig) - print*,'fm(ig,l+1)',fm(ig,l+1) - print*,'fm(ig,l)',fm(ig,l) - abort_message = 'probleme dans thermcell_flux' - CALL abort_physic (modname,abort_message,1) - endif - 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 - endif -! -! 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 -!----------------------------------------------------------------------- - - if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then - if (prt_level.ge.10) then - 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 - endif - - if (prt_level.ge.10) & - & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & - & ptimestep,masse,entr,detr,fm,'fin') - - - return - end - - subroutine printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & - & ptimestep,masse,entr,detr,fm,descr) - - implicit none - - integer ngrid,klev,lunout,igout,l,lm - - integer lmax(klev),lalim(klev) - real ptimestep,masse(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev) - real fm(ngrid,klev+1),f(ngrid) - - character*3 descr - - character (len=20) :: modname='thermcell_flux' - character (len=80) :: abort_message - - lm=lmax(igout)+5 - if(lm.gt.klev) lm=klev - - print*,'Impression jusque lm=',lm - - write(lunout,*) 'Dans thermcell_flux '//descr - write(lunout,*) 'flux base ',f(igout) - write(lunout,*) 'lmax ',lmax(igout) - write(lunout,*) 'lalim ',lalim(igout) - write(lunout,*) 'ig= ',igout - write(lunout,'(a3,4a14)') 'l','M','E','D','F' - write(lunout,'(i4,4e14.4)') (l,masse(igout,l)/ptimestep, & - & entr(igout,l),detr(igout,l) & - & ,fm(igout,l+1),l=1,lm) - - - do l=lmax(igout)+1,klev - if (abs(entr(igout,l))+abs(detr(igout,l))+abs(fm(igout,l)).gt.0.) then - print*,'cas 1 : igout,l,lmax(igout)',igout,l,lmax(igout) - print*,'entr(igout,l)',entr(igout,l) - print*,'detr(igout,l)',detr(igout,l) - print*,'fm(igout,l)',fm(igout,l) - abort_message = '' - CALL abort_physic (modname,abort_message,1) - endif - enddo - - return - end - Index: /LMDZ6/trunk/libf/phylmd/thermcell_ini_mod.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_ini_mod.F90 (revision 4090) +++ /LMDZ6/trunk/libf/phylmd/thermcell_ini_mod.F90 (revision 4090) @@ -0,0 +1,102 @@ +MODULE thermcell_ini_mod +IMPLICIT NONE + +save + integer :: dvdq=1,dqimpl=-1,prt_level=0,lunout + real RG,RD,RCPD,RKAPPA,RLVTT,RLvCp,RETV + real :: r_aspect_thermals,tau_thermals,fact_thermals_ed_dz + integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure + + +!$OMP THREADPRIVATE(dvdq,dqimpl,prt_level,lunout) +!$OMP THREADPRIVATE(RG,RD,RCPD,RKAPPA,RLVTT,RLvCp) +!$OMP THREADPRIVATE(r_aspect_thermals,tau_thermals,fact_thermals_ed_dz) +!$OMP THREADPRIVATE(iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure) + + REAL, SAVE :: fact_epsilon=0.002 + REAL, SAVE :: betalpha=0.9 + REAL, SAVE :: afact=2./3. + REAL, SAVE :: fact_shell=1. + REAL,SAVE :: detr_min=1.e-5 + REAL,SAVE :: entr_min=1.e-5 + REAL,SAVE :: detr_q_coef=0.012 + REAL,SAVE :: detr_q_power=0.5 + REAL,SAVE :: mix0=0. + INTEGER,SAVE :: thermals_flag_alim=0 + +!$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell) +!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power) +!$OMP THREADPRIVATE( mix0, thermals_flag_alim) + + +CONTAINS + +SUBROUTINE thermcell_ini(iflag_thermals,prt_level_in,tau_thermals_in,lunout_in, & + & RG_in,RD_in,RCPD_in,RKAPPA_in,RLVTT_in,RETV_in) + + USE ioipsl_getin_p_mod, ONLY : getin_p + +integer, intent(in) :: iflag_thermals,prt_level_in,lunout_in +real, intent(in) :: RG_in,RD_in,RCPD_in,RKAPPA_in,RLVTT_in,RETV_in,tau_thermals_in + +print*,'thermcell_ini' + if (iflag_thermals==15.or.iflag_thermals==16) then + dvdq=0 + dqimpl=-1 + else + dvdq=1 + dqimpl=1 + endif + prt_level=prt_level_in + RG=RG_in + RD=RD_in + RCPD=RCPD_in + RKAPPA=RKAPPA_in + RLVTT=RLVTT_in + RLvCp = RLVTT/RCPD + RETV=RETV_in + tau_thermals=tau_thermals_in + lunout=lunout_in + + +!===================================================================== +! a la fois les vieilles param et thermcell_main : +!===================================================================== + + r_aspect_thermals=2. + CALL getin_p('r_aspect_thermals',r_aspect_thermals) + + tau_thermals = 0. + CALL getin_p('tau_thermals',tau_thermals) + + fact_thermals_ed_dz = 0.1 + CALL getin_p('fact_thermals_ed_dz',fact_thermals_ed_dz) + + fact_thermals_ed_dz = 0.1 + CALL getin_p('fact_thermals_ed_dz',fact_thermals_ed_dz) + + iflag_thermals_ed = 0 + CALL getin_p('iflag_thermals_ed',iflag_thermals_ed) + + iflag_thermals_optflux = 0 + CALL getin_p('iflag_thermals_optflux',iflag_thermals_optflux) + + iflag_thermals_closure = 1 + CALL getin_p('iflag_thermals_closure',iflag_thermals_closure) + + CALL getin_p('thermals_fact_epsilon',fact_epsilon) + CALL getin_p('thermals_betalpha',betalpha) + CALL getin_p('thermals_afact',afact) + CALL getin_p('thermals_fact_shell',fact_shell) + CALL getin_p('thermals_detr_min',detr_min) + CALL getin_p('thermals_entr_min',entr_min) + CALL getin_p('thermals_detr_q_coef',detr_q_coef) + CALL getin_p('thermals_detr_q_power',detr_q_power) + CALL getin_p('thermals_mix0',mix0) + CALL getin_p('thermals_flag_alim',thermals_flag_alim) + + +!===================================================================== + +END SUBROUTINE thermcell_ini +END MODULE thermcell_ini_mod Index: DZ6/trunk/libf/phylmd/thermcell_init.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_init.F90 (revision 4089) +++ (revision ) @@ -1,59 +1,0 @@ -! -! $Header$ -! - SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlay,zlev, & - & lalim,lmin,alim_star,alim_star_tot,lev_out) - -!---------------------------------------------------------------------- -!thermcell_init: calcul du profil d alimentation du thermique -!---------------------------------------------------------------------- - USE print_control_mod, ONLY: lunout - IMPLICIT NONE -#include "thermcell.h" - - INTEGER l,ig -!arguments d entree - INTEGER ngrid,nlay - REAL ztv(ngrid,nlay) - REAL zlay(ngrid,nlay) - REAL zlev(ngrid,nlay+1) -!arguments de sortie - INTEGER lalim(ngrid) - INTEGER lmin(ngrid) - REAL alim_star(ngrid,nlay) - REAL alim_star_tot(ngrid) - integer lev_out ! niveau pour les print - - REAL zzalim(ngrid) -!CR: ponderation entrainement des couches instables -!def des alim_star tels que alim=f*alim_star - - - write(lunout,*)'THERM INIT V20C ' - - alim_star_tot(:)=0. - alim_star(:,:)=0. - lmin(:)=1 - lalim(:)=1 - - do l=1,nlay-1 - do ig=1,ngrid - if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then - alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & - & *sqrt(zlev(ig,l+1)) - lalim(:)=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. - - return - end Index: /LMDZ6/trunk/libf/phylmd/thermcell_old.h =================================================================== --- /LMDZ6/trunk/libf/phylmd/thermcell_old.h (revision 4090) +++ /LMDZ6/trunk/libf/phylmd/thermcell_old.h (revision 4090) @@ -0,0 +1,4 @@ + + real,parameter :: r_aspect_thermals=2. + real,parameter :: l_mix_thermals=30. + integer,parameter :: w2di_thermals=0 Index: /LMDZ6/trunk/libf/phylmdiso/alpale.h =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/alpale.h (revision 4090) +++ /LMDZ6/trunk/libf/phylmdiso/alpale.h (revision 4090) @@ -0,0 +1,1 @@ +link ../phylmd/alpale.h Index: DZ6/trunk/libf/phylmdiso/sisvat =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/sisvat (revision 4089) +++ (revision ) @@ -1,1 +1,0 @@ -link ../phylmd/sisvat/ Index: DZ6/trunk/libf/phylmdiso/thermcell.h =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/thermcell.h (revision 4089) +++ (revision ) @@ -1,1 +1,0 @@ -link ../phylmd/thermcell.h Index: DZ6/trunk/libf/phylmdiso/thermcellV0_main.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/thermcellV0_main.F90 (revision 4089) +++ (revision ) @@ -1,1 +1,0 @@ -link ../phylmd/thermcellV0_main.F90 Index: DZ6/trunk/libf/phylmdiso/thermcell_condens.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/thermcell_condens.F90 (revision 4089) +++ (revision ) @@ -1,1 +1,0 @@ -link ../phylmd/thermcell_condens.F90 Index: DZ6/trunk/libf/phylmdiso/thermcell_flux.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/thermcell_flux.F90 (revision 4089) +++ (revision ) @@ -1,1 +1,0 @@ -link ../phylmd/thermcell_flux.F90 Index: /LMDZ6/trunk/libf/phylmdiso/thermcell_ini_mod.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/thermcell_ini_mod.F90 (revision 4090) +++ /LMDZ6/trunk/libf/phylmdiso/thermcell_ini_mod.F90 (revision 4090) @@ -0,0 +1,1 @@ +link ../phylmd/thermcell_ini_mod.F90 Index: DZ6/trunk/libf/phylmdiso/thermcell_init.F90 =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/thermcell_init.F90 (revision 4089) +++ (revision ) @@ -1,1 +1,0 @@ -link ../phylmd/thermcell_init.F90 Index: /LMDZ6/trunk/libf/phylmdiso/thermcell_old.h =================================================================== --- /LMDZ6/trunk/libf/phylmdiso/thermcell_old.h (revision 4090) +++ /LMDZ6/trunk/libf/phylmdiso/thermcell_old.h (revision 4090) @@ -0,0 +1,1 @@ +link ../phylmd/thermcell_old.h