! ! ! SUBROUTINE thermcell_main(ngrid,nlay,nq,ptimestep,firstcall, & pplay,pplev,pphi,zpopsk, & pu,pv,pt,pq, & pduadj,pdvadj,pdtadj,pdqadj, & f0,fm0,entr0,detr0,zw2,fraca, & zqta,zqla,ztv,ztva,zhla,zhl,zqsa, & lmin,lmix,lmax) !=============================================================================== ! 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 ; impl = 1 : implicit ; impl =-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" ! ! Major changes 2018-19 (AB alexandre.boissinot@lmd.jussieu.fr) ! New detr and entre formulae (no longer alimentation) ! lmin can be greater than 1 ! Mix every tracer (EN COURS) ! Old version of thermcell_dq is removed ! Alternative version thermcell_dv2 is removed ! !=============================================================================== USE thermcell_mod USE print_control_mod, ONLY: prt_level IMPLICIT NONE !=============================================================================== ! Declaration !=============================================================================== ! Inputs: ! ------- INTEGER ngrid, nlay, nq REAL ptimestep REAL pplay(ngrid,nlay) ! Layer pressure REAL pplev(ngrid,nlay+1) ! Level pressure REAL pphi(ngrid,nlay) ! Geopotential REAL pu(ngrid,nlay) ! Zonal wind REAL pv(ngrid,nlay) ! Meridional wind REAL pt(ngrid,nlay) ! Temperature REAL pq(ngrid,nlay,nq) ! Tracers mass mixing ratio LOGICAL firstcall ! Outputs: ! -------- REAL pduadj(ngrid,nlay) ! u convective variations REAL pdvadj(ngrid,nlay) ! v convective variations REAL pdtadj(ngrid,nlay) ! t convective variations REAL pdqadj(ngrid,nlay,nq) ! q convective variations REAL f0(ngrid) ! mass flux norm (after possible time relaxation) REAL fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation) REAL entr0(ngrid,nlay) ! entrainment (after possible time relaxation) REAL detr0(ngrid,nlay) ! detrainment (after possible time relaxation) ! Local: ! ------ INTEGER ig, k, l, iq INTEGER lmax(ngrid) ! Highest layer reached by the plume INTEGER lmix(ngrid) ! Layer in which plume vertical speed is maximal INTEGER lmin(ngrid) ! First unstable layer REAL zmix(ngrid) ! Altitude of maximal vertical speed REAL zmax(ngrid) ! Maximal altitudes where plumes are active REAL zmin(ngrid) ! Minimal altitudes where plumes are active REAL zlay(ngrid,nlay) ! Layers altitudes REAL zlev(ngrid,nlay+1) ! Levels altitudes REAL rho(ngrid,nlay) ! Layers densities REAL rhobarz(ngrid,nlay) ! Levels densities REAL masse(ngrid,nlay) ! Layers masses REAL zpopsk(ngrid,nlay) ! Exner function REAL zu(ngrid,nlay) ! u environment REAL zv(ngrid,nlay) ! v environment REAL zt(ngrid,nlay) ! TR environment REAL zqt(ngrid,nlay) ! qt environment REAL zql(ngrid,nlay) ! ql environment REAL zhl(ngrid,nlay) ! TP environment REAL ztv(ngrid,nlay) ! TRPV environment REAL zqs(ngrid,nlay) ! qsat environment REAL zua(ngrid,nlay) ! u plume REAL zva(ngrid,nlay) ! v plume REAL zta(ngrid,nlay) ! TR plume REAL zqla(ngrid,nlay) ! qv plume REAL zqta(ngrid,nlay) ! qt plume REAL zhla(ngrid,nlay) ! TP plume REAL ztva(ngrid,nlay) ! TRPV plume REAL zqsa(ngrid,nlay) ! qsat plume REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt) REAL f_star(ngrid,nlay+1) ! Normalized mass flux REAL entr_star(ngrid,nlay) ! Normalized entrainment (E* = e* dz) REAL detr_star(ngrid,nlay) ! Normalized detrainment (D* = d* dz) REAL fm(ngrid,nlay+1) ! Mass flux REAL entr(ngrid,nlay) ! Entrainment (E = e dz) REAL detr(ngrid,nlay) ! Detrainment (D = d dz) REAL f(ngrid) ! Mass flux norm REAL lambda ! Time relaxation coefficent REAL fraca(ngrid,nlay+1) ! Updraft fraction REAL linter(ngrid) ! Level (continuous) of maximal vertical speed REAL wmax(ngrid) ! Maximal vertical speed REAL zw2(ngrid,nlay+1) ! Plume vertical speed REAL zdthladj(ngrid,nlay) ! Potential temperature variations REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() !=============================================================================== ! Initialization !=============================================================================== IF (firstcall) THEN fm0(:,:) = 0. entr0(:,:) = 0. detr0(:,:) = 0. ENDIF f_star(:,:) = 0. entr_star(:,:) = 0. detr_star(:,:) = 0. f(:) = 0. fm(:,:) = 0. entr(:,:) = 0. detr(:,:) = 0. lmax(:) = 1 lmix(:) = 1 lmin(:) = 1 pduadj(:,:) = 0.0 pdvadj(:,:) = 0.0 pdtadj(:,:) = 0.0 pdqadj(:,:,:) = 0.0 DO ig=1,ngrid ! AB: Careful: Hard-coded value from Earth version! ! f0(ig) = max(f0(ig), 1.e-2) ! AB: No pescribed minimal value for f0 f0(ig) = max(f0(ig), 0.) ENDDO !=============================================================================== ! Environment settings !=============================================================================== !------------------------------------------------------------------------------- ! Calcul de T,q,ql a partir de Tl et qt dans l environnement !------------------------------------------------------------------------------- CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev, & & zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs) !------------------------------------------------------------------------------- ! Levels and layers altitudes !------------------------------------------------------------------------------- DO l=2,nlay zlev(:,l) = 0.5 * (pphi(:,l) + pphi(:,l-1)) / RG ENDDO zlev(:,1) = 0. zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG DO l=1,nlay zlay(:,l) = pphi(:,l) / RG ENDDO !------------------------------------------------------------------------------- ! Levels and layers densities !------------------------------------------------------------------------------- rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:)) IF (prt_level.ge.10) THEN print *, 'WARNING: density in the first layer is equal to density at the first level!' print *, 'rhobarz(:,1)', rhobarz(:,1) print *, 'rho(:,1) ', rho(:,1) ENDIF rhobarz(:,1) = rho(:,1) DO l=2,nlay rhobarz(:,l) = 0.5 * (rho(:,l) + rho(:,l-1)) ENDDO !------------------------------------------------------------------------------- ! Layers masses !------------------------------------------------------------------------------- DO l=1,nlay masse(:,l) = (pplev(:,l) - pplev(:,l+1)) / RG ENDDO !=============================================================================== ! Explicative schemes !=============================================================================== !------------------------------------------------------------------------------- ! Thermal plume variables !------------------------------------------------------------------------------- ! top of the model ! =========================== ! ! --------------------------- ! _ ! ----- F_lmax+1=0 ------zmax \ ! lmax | ! ------F_lmax>0------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | ! ------------------wmax,zmix | ! lmix | ! --------------------------- | ! | ! --------------------------- | ! | E, D ! --------------------------- | ! | ! --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca ! zt, zu, zv, zo, rho | ! --------------------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | ! ------F_lmin+1>0----------- | ! lmin | ! ----- F_lmin=0 ------------ _/ ! ! --------------------------- ! ! =========================== ! bottom of the model !------------------------------------------------------------------------------- ! Zoom on layers k and k-1 !------------------------------------------------------------------------------- ! | /|\ | | ! |---- | F_k+1 -----------|--------------------------| level k+1 ! | | w_k+1 | | ! | --|--> D_k | ! | | | layer k ! | <--|-- E_k | ! | /|\ | | ! |---- | F_k ----------|-----------------------------| level k ! | | w_k | | ! | --|--> D_k-1 | ! | | | layer k-1 ! | <--|-- E_k-1 | ! | /|\ | | ! |---- | F_k-1 -----|--------------------------------| level k-1 ! | w_k-1 ! 0 fraca 1 ! \__________________/ \______________________________/ ! plume (fraca) environment (1-fraca) !=============================================================================== ! Thermal plumes computation !=============================================================================== !------------------------------------------------------------------------------- ! Thermal plumes speeds, fluxes, tracers and temperatures !------------------------------------------------------------------------------- CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv, & & zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, & & detr_star,entr_star,f_star, & & ztva,zhla,zqla,zqta,zta,zqsa, & & zw2,lmix,lmin) !------------------------------------------------------------------------------- ! Thermal plumes characteristics: zmax, zmix, wmax !------------------------------------------------------------------------------- ! AB: Careful, zw2 became its square root in thermcell_height! CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2, & & zlev,zmin,zmix,zmax,wmax,f_star) !=============================================================================== ! Closure and mass fluxes computation !=============================================================================== !------------------------------------------------------------------------------- ! Closure !------------------------------------------------------------------------------- CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & & lmax,entr_star,zmin,zmax,wmax,f) 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 print *, 'ERROR: mass flux norm is not positive!' print *, 'f0 =', f0(1) CALL abort ENDIF !------------------------------------------------------------------------------- ! Mass fluxes !------------------------------------------------------------------------------- CALL thermcell_flux(ngrid,nlay,ptimestep,masse, & & lmin,lmax,entr_star,detr_star, & & f,rhobarz,zlev,zw2,fm,entr,detr,zqla) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 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 sur Terre). !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 !------------------------------------------------------------------------------- ! Updraft fraction !------------------------------------------------------------------------------- DO ig=1,ngrid fraca(ig,1) = 0. fraca(ig,nlay+1) = 0. ENDDO DO l=2,nlay DO ig=1,ngrid IF (zw2(ig,l) > 0.) THEN fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l)) ELSE fraca(ig,l) = 0. ENDIF ENDDO ENDDO !=============================================================================== ! Transport vertical !=============================================================================== !------------------------------------------------------------------------------- ! Calcul du transport vertical de la temperature potentielle !------------------------------------------------------------------------------- CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & & zhl,zdthladj,dummy,lmin) DO l=1,nlay DO ig=1,ngrid pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l) ENDDO ENDDO !------------------------------------------------------------------------------- ! Calcul du transport vertical des traceurs !------------------------------------------------------------------------------- DO iq=1,nq CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & & pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin) ENDDO !------------------------------------------------------------------------------- ! Calcul du transport vertical du moment horizontal !------------------------------------------------------------------------------- IF (dvimpl) THEN CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,fraca, & & zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva) ELSE CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & & zu,pduadj,zua,lmin) CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & & zv,pdvadj,zva,lmin) ENDIF RETURN END