! ! ! SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep & & ,pplay,pplev,pphi,debut & & ,pu,pv,pt,po & & ,pduadj,pdvadj,pdtadj,pdoadj & & ,f0,fm0,entr0,detr0 & & ,zqta,zqla,ztv,ztva,ztla,zthl,zqsatth & & ,zmax0,zw2,fraca & & ,lmin,lmix,lalim,lmax & & ,zpspsk,ratqscth,ratqsdiff & & ,Ale_bl,Alp_bl,lalim_conv,wght_th & !!! nrlmd le 10/04/2012 & ,pbl_tke,pctsrf,omega,airephy & & ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 & & ,n2,s2,ale_bl_stat & & ,therm_tke_max,env_tke_max & & ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke & & ,alp_bl_conv,alp_bl_stat) !!! fin nrlmd le 10/04/2012 !============================================================================== ! 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" ! !============================================================================== USE dimphy !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! AB : why use dimphy to get klon, klev while subroutine arguments ngrid and ! nlay are the same? !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE thermcell_mod USE print_control_mod, ONLY: lunout,prt_level USE planete_mod, ONLY: preff IMPLICIT NONE !============================================================================== ! Declaration !============================================================================== ! inputs: ! ------- INTEGER itap INTEGER ngrid INTEGER nlay REAL ptimestep REAL pplay(ngrid,nlay) REAL pplev(ngrid,nlay+1) REAL pphi(ngrid,nlay) REAL pt(ngrid,nlay) ! REAL pu(ngrid,nlay) ! REAL pv(ngrid,nlay) ! REAL po(ngrid,nlay) ! LOGICAL debut ! outputs: ! -------- REAL pdtadj(ngrid,nlay) ! t convective variations REAL pduadj(ngrid,nlay) ! u convective variations REAL pdvadj(ngrid,nlay) ! v convective variations REAL pdoadj(ngrid,nlay) ! water convective variations REAL fm0(klon,klev+1) ! mass flux (after possible time relaxation) REAL entr0(klon,klev) ! entrainment (after possible time relaxation) REAL detr0(klon,klev) ! detrainment (after possible time relaxation) REAL f0(klon) ! mass flux norm (after possible time relaxation) ! local: ! ------ INTEGER, save :: dvimpl=1 !$OMP THREADPRIVATE(dvimpl) INTEGER, save :: dqimpl=-1 !$OMP THREADPRIVATE(dqimpl) 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,ierr INTEGER lmix_bis(klon) INTEGER lmax(klon) ! INTEGER lmin(klon) ! INTEGER lalim(klon) ! INTEGER lmix(klon) ! REAL linter(klon) REAL zmix(klon) REAL zmax(klon) REAL zw2(klon,klev+1) REAL zw_est(klon,klev+1) REAL zmax_sec(klon) REAL zmax0(klon) REAL zlay(klon,klev) ! layers altitude REAL zlev(klon,klev+1) ! levels altitude REAL rho(klon,klev) ! layers density REAL rhobarz(klon,klev) ! levels density REAL deltaz(klon,klev) ! layers heigth REAL masse(klon,klev) ! layers mass REAL zpspsk(klon,klev) ! Exner function REAL zu(klon,klev) ! environment zonal wind REAL zv(klon,klev) ! environment meridional wind REAL zo(klon,klev) ! environment water tracer REAL zva(klon,klev) ! plume zonal wind REAL zua(klon,klev) ! plume meridional wind REAL zoa(klon,klev) ! plume water tracer REAL zt(klon,klev) ! T environment REAL zh(klon,klev) ! T,TP environment REAL zthl(klon,klev) ! TP environment REAL ztv(klon,klev) ! TPV environment ? REAL zl(klon,klev) ! ql environment REAL zta(klon,klev) ! REAL zha(klon,klev) ! TRPV plume REAL ztla(klon,klev) ! TP plume REAL ztva(klon,klev) ! TRPV plume REAL ztva_est(klon,klev) ! TRPV plume (temporary) REAL zqla(klon,klev) ! qv plume REAL zqta(klon,klev) ! qt plume REAL wmax(klon) ! maximal vertical speed REAL wmax_tmp(klon) ! temporary maximal vertical speed REAL wmax_sec(klon) ! maximal vertical speed if dry case REAL fraca(klon,klev+1) ! updraft fraction REAL f_star(klon,klev+1) ! normalized mass flux REAL entr_star(klon,klev) ! normalized entrainment REAL detr_star(klon,klev) ! normalized detrainment REAL alim_star_tot(klon) ! integrated alimentation REAL alim_star(klon,klev) ! normalized alimentation REAL alim_star_clos(klon,klev) ! closure alimentation REAL fm(klon,klev+1) ! mass flux REAL entr(klon,klev) ! entrainment REAL detr(klon,klev) ! detrainment REAL f(klon) ! mass flux norm REAL zdthladj(klon,klev) ! REAL lambda ! time relaxation coefficent REAL zsortie(klon,klev) REAL zsortie1d(klon) REAL susqr2pi, Reuler REAL zf REAL zf2 REAL thetath2(klon,klev) REAL wth2(klon,klev) REAL wth3(klon,klev) REAL q2(klon,klev) ! FH : probleme de dimensionnement avec l'allocation dynamique ! common/comtherm/thetath2,wth2 real wq(klon,klev) real wthl(klon,klev) real wthv(klon,klev) real ratqscth(klon,klev) real var real vardiff real ratqsdiff(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) ! FH/IM : save f0 real zlevinter(klon) real seuil !!! nrlmd le 10/04/2012 !------Entrees real pbl_tke(klon,klev+1,nbsrf) real pctsrf(klon,nbsrf) real omega(klon,klev) real airephy(klon) !------Sorties real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon) real therm_tke_max0(klon),env_tke_max0(klon) real n2(klon),s2(klon) real ale_bl_stat(klon) real therm_tke_max(klon,klev),env_tke_max(klon,klev) real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon) !------Local integer nsrf real rhobarz0(klon) ! Densite au LCL logical ok_lcl(klon) ! Existence du LCL des thermiques integer klcl(klon) ! Niveau du LCL real interp(klon) ! Coef d'interpolation pour le LCL !------Triggering real Su ! Surface unite: celle d'un updraft elementaire parameter(Su=4e4) real hcoef ! Coefficient directeur pour le calcul de s2 parameter(hcoef=1) real hmincoef ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 parameter(hmincoef=0.3) real eps1 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd) parameter(eps1=0.3) real hmin(ngrid) ! Ordonnee a l'origine pour le calcul de s2 real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl) real zmax_moy_coef parameter(zmax_moy_coef=0.33) real depth(klon) ! Epaisseur moyenne du cumulus real w_max(klon) ! Vitesse max statistique real s_max(klon) !------Closure real pbl_tke_max(klon,klev) ! Profil de TKE moyenne real pbl_tke_max0(klon) ! TKE moyenne au LCL real w_ls(klon,klev) ! Vitesse verticale grande echelle (m/s) real coef_m ! On considere un rendement pour alp_bl_fluct_m parameter(coef_m=1.) real coef_tke ! On considere un rendement pour alp_bl_fluct_tke parameter(coef_tke=1.) !!! fin nrlmd le 10/04/2012 ! Nouvelles variables pour la convection real Ale_bl(klon) real Alp_bl(klon) real alp_int(klon),dp_int(klon),zdp real ale_int(klon) integer n_int(klon) real fm_tot(klon) real wght_th(klon,klev) integer lalim_conv(klon) CHARACTER*2 str2 CHARACTER*10 str10 CHARACTER (len=20) :: modname='thermcell_main' CHARACTER (len=80) :: abort_message EXTERNAL SCOPY !============================================================================== ! Initialization !============================================================================== seuil = 0.25 IF (debut) THEN IF (iflag_thermals==15.or.iflag_thermals==16) THEN dvimpl = 0 dqimpl = -1 ELSE dvimpl = 1 dqimpl = 1 ENDIF fm0(:,:) = 0. entr0(:,:) = 0. detr0(:,:) = 0. ENDIF fm(:,:) = 0. entr(:,:) = 0. detr(:,:) = 0. IF (ngrid.ne.klon) THEN print *, 'ERROR: ngrid and klon are different!' print *, 'ngrid =', ngrid print *, 'klon =', klon ENDIF DO ig=1,klon f0(ig) = max(f0(ig), 1.e-2) zmax0(ig) = max(zmax0(ig),40.) ENDDO IF (prt_level.ge.20) then DO ig=1,ngrid print *, 'ig,f0', ig, f0(ig) ENDDO ENDIF wmax_tmp(:) = 0. !------------------------------------------------------------------------------ ! 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) !------------------------------------------------------------------------------ ! ! -------------------- ! ! ! + + + + + + + + + + + ! ! ! 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 !------------------------------------------------------------------------------ ! Calcul des densites !------------------------------------------------------------------------------ rho(:,:) = pplay(:,:) / (zpspsk(:,:) * RD * ztv(:,:)) IF (prt_level.ge.10) THEN write(lunout,*) 'WARNING: thermcell_main rhobarz(:,1)=rho(:,1)' ENDIF 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 | | ! ----- 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 ! 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. !------------------------------------------------------------------------------ entr_star(:,:) = 0. detr_star(:,:) = 0. alim_star(:,:) = 0. alim_star_tot(:) = 0. lmin(:) = 1 !============================================================================== ! Calcul du melange et des variables dans le thermique !============================================================================== CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl, & & po,zl,rhobarz,zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & & lalim,f0,detr_star,entr_star,f_star,ztva, & & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis, & & lmin,lev_out,lunout1,igout) CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ') !============================================================================== ! Calcul des caracteristiques du thermique:zmax,zmix,wmax !============================================================================== CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, & & zlev,lmax,zmax,zmax0,zmix,wmax,f_star, & & lev_out) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! AB : WARNING: zw2 became its square root in thermcell_height! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ') !============================================================================== ! Closure and mass flux determining !============================================================================== CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & & lalim,lmin,zmax_sec,wmax_sec,lev_out) CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ') !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 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(:,:) !------------------------------------------------------------------------------ ! Closure (dry if iflag_thermals_closure=1, moist if iflag_thermals_closure=2) !------------------------------------------------------------------------------ IF (iflag_thermals_closure.eq.1) THEN CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & & lalim,alim_star_clos,f_star, & & zmax_sec,wmax_sec,f,lev_out) ELSEIF (iflag_thermals_closure.eq.2) THEN CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & & lalim,alim_star,f_star, & & zmax,wmax,f,lev_out) ELSE print *, 'ERROR: no closure had been selected!' call abort ENDIF 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.)' write(*,*) 'f0 =', f0(1) CALL abort_physic(modname,abort_message,1) ENDIF !============================================================================== ! Deduction des flux !============================================================================== CALL thermcell_flux(ngrid,nlay,ptimestep,masse, & & lalim,lmin,lmax,alim_star,entr_star,detr_star, & & f,rhobarz,zlev,zw2,fm,entr, & & detr,zqla,lev_out,lunout1,igout) CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') CALL test_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 !============================================================================== ! Calcul du transport vertical (de qt et tp) !============================================================================== CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & zthl,zdthladj,zta,lmin,lev_out) CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & po,pdoadj,zoa,lmin,lev_out) DO l=1,nlay DO ig=1,ngrid pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l) ENDDO ENDDO !============================================================================== ! 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 !============================================================================== IF (dvimpl.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) ELSE !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Calcul purement conservatif pour le transport de V !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & zu,pduadj,zua,lmin,lev_out) CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & zv,pdvadj,zva,lmin,lev_out) ENDIF !============================================================================== ! Calculs de diagnostiques pour les sorties !============================================================================== IF (sorties) THEN !------------------------------------------------------------------------------ ! Calcul du niveau de condensation !------------------------------------------------------------------------------ if (prt_level.ge.1) print*,'14a OK convect8' 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 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 write(*,*) 'ERROR: thermal plumes rise the model top' CALL abort 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 !------------------------------------------------------------------------------ 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' 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.-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 !------------------------------------------------------------------------------ ! 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.-po(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 CALL thermcell_alp(ngrid,nlay,ptimestep, & & pplay,pplev, & & fm0,entr0,lmax, & & Ale_bl,Alp_bl,lalim_conv,wght_th, & & zw2,fraca, & & pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax, & & pbl_tke,pctsrf,omega,airephy, & & zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0, & & n2,s2,ale_bl_stat, & & therm_tke_max,env_tke_max, & & alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke, & & alp_bl_conv,alp_bl_stat) !------------------------------------------------------------------------------ ! Calcul du ratqscdiff !------------------------------------------------------------------------------ var = 0. vardiff = 0. ratqsdiff(:,:) = 0. DO l=1,klev 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,klev 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) / (po(ig,l) * 1000.) ENDDO ENDDO ENDIF RETURN END !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po, & ztva,zqla,f_star,zw2,comment) USE print_control_mod, ONLY: prt_level IMPLICIT NONE !============================================================================== ! Declaration !============================================================================== ! inputs: ! ------- INTEGER klon INTEGER klev INTEGER long(klon) REAL pplev(klon,klev+1) REAL 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) REAL seuil CHARACTER*21 comment ! local: ! ------ INTEGER i, k !============================================================================== ! Test !============================================================================== IF (prt_level.ge.1) THEN write(*,*) 'WARNING: in 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 thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, & & rg,pplev,therm_tke_max) USE print_control_mod, 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 ! ! Transport de la TKE par le thermique moyen pour la fermeture en ALP ! On transporte pbl_tke pour donner therm_tke !============================================================================== integer ngrid,nlay,nsrf real ptimestep real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1) real entr0(ngrid,nlay),rg real therm_tke_max(ngrid,nlay) real detr0(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) real zzm integer ig,k integer isrf 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. !!! nrlmd le 16/09/2010 ! calcul de la valeur dans les ascendances ! do ig=1,ngrid ! qa(ig,1)=q(ig,1) ! enddo !!! !do isrf=1,nsrf ! q(:,:)=therm_tke(:,:,isrf) 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 END