! $Id: lmdz_thermcell_plume_6A.f90 5116 2024-07-24 12:54:37Z abarral $ MODULE lmdz_thermcell_plume_6A CONTAINS SUBROUTINE thermcell_plume_6A(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) !-------------------------------------------------------------------------- !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance !-------------------------------------------------------------------------- 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 integer, intent(out), dimension(ngrid) :: lmix_bis real, intent(out), dimension(ngrid) :: linter REAL zdw2, zdw2bis REAL zw2modif REAL zw2fact, zw2factbis REAL, dimension(ngrid, nlay) :: zeps REAL, dimension(ngrid) :: wmaxa INTEGER ig, l, k, lt, it, lm integer nbpb real, dimension(ngrid, nlay) :: detr real, dimension(ngrid, nlay) :: entr real, dimension(ngrid, nlay + 1) :: wa_moy real, dimension(ngrid, nlay) :: ztv_est real, dimension(ngrid) :: ztemp, zqsat real, dimension(ngrid, nlay) :: zqla_est real, dimension(ngrid, nlay) :: zta_est real, dimension(ngrid, nlay) :: zbuoy, gamma, zdqt real zdz, zalpha, zw2m real, dimension(ngrid, nlay) :: zbuoyjam, zdqtjam real zbuoybis, 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 coefc 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. ! Initialisation a 0 en cas de sortie dans replay zqsat(:) = 0. zta_est(:, :) = 0. zdqt(:, :) = 0. zdqtjam(:, :) = 0. c2(:, :) = 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 !========================================================================= !-------------------------------------------------- if (iflag_thermals_ed<8) THEN !-------------------------------------------------- !AJ052014: J'ai remplac?? la boucle do par un do while ! afin de faire moins de calcul dans la boucle !-------------------------------------------------- do while (zlmelup>zltup) lt = lt + 1 zlt = zlev(ig, lt) zdz3 = zlev(ig, lt + 1) - zlt zltdwn = zlt - zdz3 / 2 zltup = zlt + zdz3 / 2 enddo !-------------------------------------------------- !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors ! on cherche o?? se trouve l'altitude d'inversion ! en calculant ztv1 (interpolation de la valeur de ! theta au niveau lt en utilisant les niveaux lt-1 et ! lt-2) et ztv2 (interpolation avec les niveaux lt+1 ! et lt+2). Si theta r??ellement calcul??e au niveau lt ! comprise entre ztv1 et ztv2, alors il y a inversion ! et on calcule son altitude zinv en supposant que ztv(lt) ! est une combinaison lineaire de ztv1 et ztv2. ! Ensuite, on calcule la flottabilite en comparant ! la temperature de la couche l a celle de l'air situe ! l+lmel plus haut, ce qui necessite de savoir quel fraction ! de cet air est au-dessus ou en-dessous de l'inversion !-------------------------------------------------- atv1 = (ztv(ig, lt - 1) - ztv(ig, lt - 2)) / (zlev(ig, lt - 1) - zlev(ig, lt - 2)) btv1 = (ztv(ig, lt - 2) * zlev(ig, lt - 1) - ztv(ig, lt - 1) * zlev(ig, lt - 2)) & / (zlev(ig, lt - 1) - zlev(ig, lt - 2)) atv2 = (ztv(ig, lt + 2) - ztv(ig, lt + 1)) / (zlev(ig, lt + 2) - zlev(ig, lt + 1)) btv2 = (ztv(ig, lt + 1) * zlev(ig, lt + 2) - ztv(ig, lt + 2) * zlev(ig, lt + 1)) & / (zlev(ig, lt + 2) - zlev(ig, lt + 1)) ztv1 = atv1 * zlt + btv1 ztv2 = atv2 * zlt + btv2 if (ztv(ig, lt)>ztv1.and.ztv(ig, lt)=zinv) THEN ztv_est(ig, l) = atv2 * zlmel + btv2 zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) & + (1. - fact_shell) * zbuoy(ig, l) elseif (zlmelup>=zinv) THEN ztv_est2 = atv2 * 0.5 * (zlmelup + zinv) + btv2 ztv_est1 = atv1 * 0.5 * (zinv + zlmeldwn) + btv1 ztv_est(ig, l) = ((zlmelup - zinv) / zdz) * ztv_est2 + ((zinv - zlmeldwn) / zdz) * ztv_est1 zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zinv) / zdz) * (ztva_est(ig, l) - & ztv_est2) / ztv_est2 + ((zinv - zlmeldwn) / zdz) * (ztva_est(ig, l) - & ztv_est1) / ztv_est1) + (1. - fact_shell) * zbuoy(ig, l) else ztv_est(ig, l) = atv1 * zlmel + btv1 zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) & + (1. - fact_shell) * zbuoy(ig, l) endif else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) THEN if (zlmeldwn>zltdwn) THEN zbuoyjam(ig, l) = fact_shell * RG * ((ztva_est(ig, l) - & ztv(ig, lt)) / ztv(ig, lt)) + (1. - fact_shell) * zbuoy(ig, l) else zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zltdwn) / zdz) * (ztva_est(ig, l) - & ztv(ig, lt)) / ztv(ig, lt) + ((zltdwn - zlmeldwn) / zdz) * (ztva_est(ig, l) - & ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + (1. - fact_shell) * zbuoy(ig, l) endif ! zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & ! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- & ! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l) ! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & ! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- & ! & po(ig,lt-1))/po(ig,lt-1)) endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) THEN else ! if (iflag_thermals_ed.lt.8) THEN 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) endif ! if (iflag_thermals_ed.lt.8) THEN !------------------------------------------------ !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 ! zdw2bis=0.5*(zdw2+zdw2bis) lm = Max(1, l - 2) ! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) & ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) ! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) ! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) ! w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* & ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & ! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) ! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact)) !-------------------------------------------------- !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l) !-------------------------------------------------- if (iflag_thermals_ed==8) THEN ! Ancienne version ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)) w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2) ! Nouvelle version Arnaud else ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)) w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2) ! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* & ! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* & ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis)) ! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact)) ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & ! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) ! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) endif if (iflag_thermals_ed<6) THEN zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l) ! fact_epsilon=0.0005/(zalpha+0.025)**0.5 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) fact_epsilon = 0.0002 / (zalpha + 0.1) 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 ! w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2) endif !-------------------------------------------------- !AJ052014: J'ai comment? ce if plus n?cessaire puisqu' !on fait max(0.0001,.....) !-------------------------------------------------- ! if (w_est(ig,l+1).lt.0.) THEN ! w_est(ig,l+1)=zw2(ig,l) ! w_est(ig,l+1)=0.0001 ! endif 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) ! zw2m=zw2(ig,l) zdz = zlev(ig, l + 1) - zlev(ig, l) zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l) ! zbuoybis=zbuoy(ig,l)+RG*0.1/300. zbuoybis = zbuoy(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) ! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & ! & afact*zbuoybis/zw2m - fact_epsilon ) ! entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha* & ! & afact*zbuoybis/zw2m - fact_epsilon ) ! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) !========================================================================= ! 4. Calcul de l'entrainement et du detrainement !========================================================================= ! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & ! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon ) ! entrbis=entr_star(ig,l) if (iflag_thermals_ed<6) THEN fact_epsilon = 0.0002 / (zalpha + 0.1) endif 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)) ! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ & ! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1) zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l) 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)) ! entr_star(ig,l)=f_star(ig,l)*zdz* ( & ! & mix0 * 0.1 / (zalpha+0.001) & ! & + MAX(entr_min, & ! & zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon + & ! & detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) ! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ & ! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1) ! entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha* & ! & afact*zbuoy(ig,l)/zw2m & ! & - 1.*fact_epsilon) ! En dessous de lalim, on prend le max de alim_star et entr_star pour ! alim_star et 0 sinon if (lalim_star(ig,l-1)) THEN ! alim_star(ig,l)=entrbis ! endif ! PRINT*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l) ! Calcul du flux montant normalise f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l) & - detr_star(ig, l) endif enddo !============================================================================ ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique !=========================================================================== activetmp(:) = active(:) .and. f_star(:, l + 1)>1.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) !!!!!!! fact_epsilon=0.002 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) ! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) & ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) ! lm=Max(1,l-2) ! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) ! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & ! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) ! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) if (iflag_thermals_ed==8) THEN zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2) else zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2) endif ! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* & ! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* & ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis)) if (iflag_thermals_ed<6) THEN zalpha = f0(ig) * f_star(ig, l) / sqrt(zw2(ig, l + 1)) / rhobarz(ig, l) ! fact_epsilon=0.0005/(zalpha+0.025)**0.5 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) fact_epsilon = 0.0002 / (zalpha + 0.1)**1 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,(zdz/zdzbis)*(exp(-zw2fact)* & ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) ! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2) endif 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 SUBROUTINE thermcell_plume_5B(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) !-------------------------------------------------------------------------- !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance ! Version conforme a l'article de Rio et al. 2010. ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin !-------------------------------------------------------------------------- USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG USE lmdz_thermcell_qsat, ONLY: thermcell_qsat IMPLICIT NONE INTEGER itap INTEGER lunout1, igout INTEGER ngrid, nlay REAL ptimestep REAL ztv(ngrid, nlay) REAL zthl(ngrid, nlay) REAL, INTENT(IN) :: po(ngrid, nlay) REAL zl(ngrid, nlay) REAL rhobarz(ngrid, nlay) REAL zlev(ngrid, nlay + 1) REAL pplev(ngrid, nlay + 1) REAL pphi(ngrid, nlay) REAL zpspsk(ngrid, nlay) REAL alim_star(ngrid, nlay) REAL f0(ngrid) INTEGER lalim(ngrid) integer lev_out ! niveau pour les print integer nbpb real alim_star_tot(ngrid) REAL ztva(ngrid, nlay) REAL ztla(ngrid, nlay) REAL zqla(ngrid, nlay) REAL zqta(ngrid, nlay) REAL zha(ngrid, nlay) REAL detr_star(ngrid, nlay) REAL coefc REAL entr_star(ngrid, nlay) REAL detr(ngrid, nlay) REAL entr(ngrid, nlay) REAL csc(ngrid, nlay) REAL zw2(ngrid, nlay + 1) REAL w_est(ngrid, nlay + 1) REAL f_star(ngrid, nlay + 1) REAL wa_moy(ngrid, nlay + 1) REAL ztva_est(ngrid, nlay) REAL zqla_est(ngrid, nlay) REAL zqsatth(ngrid, nlay) REAL zta_est(ngrid, nlay) REAL zbuoyjam(ngrid, nlay) REAL ztemp(ngrid), zqsat(ngrid) REAL zdw2 REAL zw2modif REAL zw2fact REAL zeps(ngrid, nlay) REAL linter(ngrid) INTEGER lmix(ngrid) INTEGER lmix_bis(ngrid) REAL wmaxa(ngrid) INTEGER ig, l, k real zdz, zbuoy(ngrid, nlay), zalpha, gamma(ngrid, nlay), zdqt(ngrid, nlay), zw2m real zbuoybis real zcor, zdelta, zcvm5, qlbef, zdz2 real betalpha, zbetalpha real eps, afact logical Zsat LOGICAL active(ngrid), activetmp(ngrid) REAL fact_gamma, fact_epsilon, fact_gamma2, fact_epsilon2 REAL c2(ngrid, nlay) Zsat = .FALSE. ! Initialisation fact_epsilon = 0.002 betalpha = 0.9 afact = 2. / 3. zbetalpha = betalpha / (1. + betalpha) ! Initialisations des variables reeles if (1==1) THEN ztva(:, :) = ztv(:, :) ztva_est(:, :) = ztva(:, :) ztla(:, :) = zthl(:, :) zqta(:, :) = po(:, :) zha(:, :) = ztva(:, :) else ztva(:, :) = 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. lalim(:) = 1 !------------------------------------------------------------------------- ! On ne considere comme actif que les colonnes dont les deux premieres ! couches sont instables. !------------------------------------------------------------------------- active(:) = ztv(:, 1)>ztv(:, 2) !------------------------------------------------------------------------- ! Definition de l'alimentation !------------------------------------------------------------------------- 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(ig) = 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. !------------------------------------------------------------------------------ ! 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)) !------------------------------------------------ !AJAM:nouveau calcul de w? !------------------------------------------------ zdz = zlev(ig, l + 1) - zlev(ig, l) zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l) zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha) zdw2 = (afact) * zbuoy(ig, l) / (fact_epsilon) w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2) if (w_est(ig, l + 1)<0.) THEN w_est(ig, l + 1) = zw2(ig, l) endif 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) ! zbuoybis=zbuoy(ig,l)+RG*0.1/300. zbuoybis = zbuoy(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) entr_star(ig, l) = f_star(ig, l) * zdz * zbetalpha * MAX(0., & afact * zbuoybis / zw2m - fact_epsilon) detr_star(ig, l) = f_star(ig, l) * zdz & * MAX(1.e-3, -afact * zbetalpha * zbuoy(ig, l) / zw2m & + 0.012 * (zdqt(ig, l) / zw2m)**0.5) ! 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) zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz) zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha) zdw2 = afact * zbuoy(ig, l) / (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 !--------------------------------------------------------------------------- !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. 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)) ! PRINT*,"linter plume", linter(ig) zw2(ig, l + 1) = 0. 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 end END MODULE lmdz_thermcell_plume_6A