MODULE lmdz_thermcell_plume ! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $ CONTAINS SUBROUTINE 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, csc, ztva, & ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter & , lev_out, lunout1, igout) ! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam) !-------------------------------------------------------------------------- ! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance ! This versions starts from a cleaning of thermcell_plume_6A (2019/01/20) ! thermcell_plume_6A is activate for flag_thermas_ed < 10 ! thermcell_plume_5B for flag_thermas_ed < 20 ! thermcell_plume for flag_thermals_ed>= 20 ! Various options are controled by the flag_thermals_ed parameter ! = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8 ! = 21 : the Jam strato-cumulus modif is not activated in detrainment ! = 29 : an other way to compute the modified buoyancy (to be tested) !-------------------------------------------------------------------------- USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim USE lmdz_thermcell_alim, ONLY: thermcell_alim USE lmdz_thermcell_qsat, ONLY: thermcell_qsat IMPLICIT NONE INTEGER, INTENT(IN) :: itap, lev_out, lunout1, igout, ngrid, nlay REAL, INTENT(IN) :: ptimestep REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: ztv REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zthl REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: po REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zl REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: rhobarz REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: zlev REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: pplev REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: pphi REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zpspsk REAL, INTENT(IN), DIMENSION(ngrid) :: f0 INTEGER, INTENT(OUT) :: lalim(ngrid) REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: alim_star REAL, INTENT(OUT), DIMENSION(ngrid) :: alim_star_tot REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: detr_star REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: entr_star REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: f_star REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: csc REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztva REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztla REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqla REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqta REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zha REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: zw2 REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: w_est REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: ztva_est REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: zqsatth INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmix(ngrid) INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmix_bis(ngrid) REAL, INTENT(OUT), DIMENSION(ngrid) :: linter(ngrid) REAL, DIMENSION(ngrid, nlay + 1) :: wa_moy REAL, DIMENSION(ngrid, nlay) :: entr, detr REAL, DIMENSION(ngrid, nlay) :: ztv_est REAL, DIMENSION(ngrid, nlay) :: zqla_est REAL, DIMENSION(ngrid, nlay) :: zta_est REAL, DIMENSION(ngrid) :: ztemp, zqsat REAL zdw2, zdw2bis REAL zw2modif REAL zw2fact, zw2factbis REAL, DIMENSION(ngrid, nlay) :: zeps REAL, DIMENSION(ngrid) :: wmaxa INTEGER ig, l, k, lt, it, lm, nbpb REAL, DIMENSION(ngrid, nlay) :: zbuoy, gamma, zdqt REAL zdz, zalpha, zw2m REAL, DIMENSION(ngrid, nlay) :: zbuoyjam, zdqtjam REAL zdz2, zdz3, lmel, entrbis, zdzbis REAL, DIMENSION(ngrid) :: d_temp REAL ztv1, ztv2, factinv, zinv, zlmel REAL zlmelup, zlmeldwn, zlt, zltdwn, zltup REAL atv1, atv2, btv1, btv2 REAL ztv_est1, ztv_est2 REAL zcor, zdelta, zcvm5, qlbef REAL zbetalpha, coefzlmel REAL eps LOGICAL Zsat LOGICAL, DIMENSION(ngrid) :: active, activetmp REAL fact_gamma, fact_gamma2, fact_epsilon2 REAL, DIMENSION(ngrid, nlay) :: c2 IF (ngrid==1) PRINT*, 'THERMCELL PLUME MODIFIE 2014/07/11' Zsat = .FALSE. ! Initialisation zbetalpha = betalpha / (1. + betalpha) ! Initialisations des variables r?elles IF (1==1) THEN ztva(:, :) = ztv(:, :) ztva_est(:, :) = ztva(:, :) ztv_est(:, :) = ztv(:, :) ztla(:, :) = zthl(:, :) zqta(:, :) = po(:, :) zqla(:, :) = 0. zha(:, :) = ztva(:, :) else ztva(:, :) = 0. ztv_est(:, :) = 0. ztva_est(:, :) = 0. ztla(:, :) = 0. zqta(:, :) = 0. zha(:, :) = 0. endif zqla_est(:, :) = 0. zqsatth(:, :) = 0. zqla(:, :) = 0. detr_star(:, :) = 0. entr_star(:, :) = 0. alim_star(:, :) = 0. alim_star_tot(:) = 0. csc(:, :) = 0. detr(:, :) = 0. entr(:, :) = 0. zw2(:, :) = 0. zbuoy(:, :) = 0. zbuoyjam(:, :) = 0. gamma(:, :) = 0. zeps(:, :) = 0. w_est(:, :) = 0. f_star(:, :) = 0. wa_moy(:, :) = 0. linter(:) = 1. ! linter(:)=1. ! Initialisation des variables entieres lmix(:) = 1 lmix_bis(:) = 2 wmaxa(:) = 0. !------------------------------------------------------------------------- ! On ne considere comme actif que les colonnes dont les deux premieres ! couches sont instables. !------------------------------------------------------------------------- active(:) = ztv(:, 1)>ztv(:, 2) d_temp(:) = 0. ! Pour activer un contraste de temperature a la base ! du panache ! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim) !------------------------------------------------------------------------------ ! Calcul dans la premiere couche ! On decide dans cette version que le thermique n'est actif que si la premiere ! couche est instable. ! Pourrait etre change si on veut que le thermiques puisse se d??clencher ! dans une couche l>1 !------------------------------------------------------------------------------ DO ig = 1, ngrid ! Le panache va prendre au debut les caracteristiques de l'air contenu ! dans cette couche. IF (active(ig)) THEN ztla(ig, 1) = zthl(ig, 1) zqta(ig, 1) = po(ig, 1) zqla(ig, 1) = zl(ig, 1) !cr: attention, prise en compte de f*(1)=1 f_star(ig, 2) = alim_star(ig, 1) zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2) & & * (zlev(ig, 2) - zlev(ig, 1)) & & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1)) w_est(ig, 2) = zw2(ig, 2) endif enddo !============================================================================== !boucle de calcul de la vitesse verticale dans le thermique !============================================================================== DO l = 2, nlay - 1 !============================================================================== ! On decide si le thermique est encore actif ou non ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test DO ig = 1, ngrid active(ig) = active(ig) & & .AND. zw2(ig, l)>1.e-10 & & .AND. f_star(ig, l) + alim_star(ig, l)>1.e-10 enddo !--------------------------------------------------------------------------- ! calcul des proprietes thermodynamiques et de la vitesse de la couche l ! sans tenir compte du detrainement et de l'entrainement dans cette ! couche ! C'est a dire qu'on suppose ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1) ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer ! avant) a l'alimentation pour avoir un calcul plus propre !--------------------------------------------------------------------------- ztemp(:) = zpspsk(:, l) * ztla(:, l - 1) CALL thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:)) DO ig = 1, ngrid ! PRINT*,'active',active(ig),ig,l IF(active(ig)) THEN zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig)) ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l) zta_est(ig, l) = ztva_est(ig, l) ztva_est(ig, l) = ztva_est(ig, l) / zpspsk(ig, l) ztva_est(ig, l) = ztva_est(ig, l) * (1. + RETV * (zqta(ig, l - 1) & - zqla_est(ig, l)) - zqla_est(ig, l)) !Modif AJAM zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l) zdz = zlev(ig, l + 1) - zlev(ig, l) lmel = fact_thermals_ed_dz * zlev(ig, l) ! lmel=0.09*zlev(ig,l) zlmel = zlev(ig, l) + lmel zlmelup = zlmel + (zdz / 2) zlmeldwn = zlmel - (zdz / 2) lt = l + 1 zlt = zlev(ig, lt) zdz3 = zlev(ig, lt + 1) - zlt zltdwn = zlt - zdz3 / 2 zltup = zlt + zdz3 / 2 !========================================================================= ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus !========================================================================= !-------------------------------------------------- lt = l + 1 zlt = zlev(ig, lt) zdz2 = zlev(ig, lt) - zlev(ig, l) DO while (lmel>zdz2) lt = lt + 1 zlt = zlev(ig, lt) zdz2 = zlt - zlev(ig, l) enddo zdz3 = zlev(ig, lt + 1) - zlt zltdwn = zlev(ig, lt) - zdz3 / 2 zlmelup = zlmel + (zdz / 2) coefzlmel = Min(1., (zlmelup - zltdwn) / zdz) zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - & ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - & ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l) !------------------------------------------------ !AJAM:nouveau calcul de w? !------------------------------------------------ zdz = zlev(ig, l + 1) - zlev(ig, l) zdzbis = zlev(ig, l) - zlev(ig, l - 1) zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l) zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha) zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha) zdw2 = afact * zbuoy(ig, l) / fact_epsilon zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon lm = Max(1, l - 2) w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2) endif enddo !------------------------------------------------- !calcul des taux d'entrainement et de detrainement !------------------------------------------------- DO ig = 1, ngrid IF (active(ig)) THEN ! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) zw2m = w_est(ig, l + 1) zdz = zlev(ig, l + 1) - zlev(ig, l) zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l) zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l) zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l) !========================================================================= ! 4. Calcul de l'entrainement et du detrainement !========================================================================= detr_star(ig, l) = f_star(ig, l) * zdz & * (mix0 * 0.1 / (zalpha + 0.001) & + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m & + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power)) IF (iflag_thermals_ed == 20) THEN entr_star(ig, l) = f_star(ig, l) * zdz * (& mix0 * 0.1 / (zalpha + 0.001) & + zbetalpha * MAX(entr_min, & afact * zbuoyjam(ig, l) / zw2m - fact_epsilon)) else entr_star(ig, l) = f_star(ig, l) * zdz * (& mix0 * 0.1 / (zalpha + 0.001) & + zbetalpha * MAX(entr_min, & afact * zbuoy(ig, l) / zw2m - fact_epsilon)) endif ! En dessous de lalim, on prend le max de alim_star et entr_star pour ! alim_star et 0 sinon IF (l1.e-10 DO ig = 1, ngrid IF (activetmp(ig)) THEN 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)) endif enddo ztemp(:) = zpspsk(:, l) * ztla(:, l) CALL thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l)) DO ig = 1, ngrid IF (activetmp(ig)) THEN ! on ecrit de maniere conservative (sat ou non) ! T = Tl +Lv/Cp ql zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l)) 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)) zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l) zdz = zlev(ig, l + 1) - zlev(ig, l) zdzbis = zlev(ig, l) - zlev(ig, l - 1) zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz) zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha) zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha) zdw2 = afact * zbuoy(ig, l) / (fact_epsilon) zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon) zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2) endif enddo IF (prt_level>=20) PRINT*, 'coucou calcul detr 460: ig, l', ig, l !=========================================================================== ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max !=========================================================================== nbpb = 0 DO ig = 1, ngrid IF (zw2(ig, l + 1)>0. .AND. zw2(ig, l + 1)<1.e-10) THEN ! stop 'On tombe sur le cas particulier de thermcell_dry' ! PRINT*,'On tombe sur le cas particulier de thermcell_plume' nbpb = nbpb + 1 zw2(ig, l + 1) = 0. linter(ig) = l + 1 endif IF (zw2(ig, l + 1)<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. !+CR:04/05/12:correction calcul linter pour calcul de zmax continu elseif (f_star(ig, l + 1)<0.) THEN linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l)) & - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l)) zw2(ig, l + 1) = 0. !fin CR:04/05/12 endif wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1)) IF (wa_moy(ig, l + 1)>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)<1.e-10) THEN lmix_bis(ig) = l + 1 endif lmix(ig) = l + 1 wmaxa(ig) = wa_moy(ig, l + 1) endif enddo IF (nbpb>0) THEN PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume' endif !========================================================================= ! FIN DE LA BOUCLE VERTICALE enddo !========================================================================= !on recalcule alim_star_tot DO ig = 1, ngrid alim_star_tot(ig) = 0. enddo DO ig = 1, ngrid DO l = 1, lalim(ig) - 1 alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l) enddo enddo IF (prt_level>=20) PRINT*, 'coucou calcul detr 470: ig, l', ig, l RETURN END END MODULE lmdz_thermcell_plume