MODULE lmdz_thermcell_alp ! $Id: thermcell_main.F90 2351 2015-08-25 15:14:59Z emillour $ CONTAINS SUBROUTINE thermcell_alp(ngrid, nlay, ptimestep & ! in , pplay, pplev & ! in , fm0, entr0, lmax & ! in , pbl_tke, pctsrf, omega, airephy & ! in , zw2, fraca & ! in , pcon, rhobarz, wth3, wmax_sec, lalim, fm, alim_star, zmax & ! in , zcong, ale_bl, alp_bl, lalim_conv, wght_th & ! out , zlcl, fraca0, w0, w_conv, therm_tke_max0, env_tke_max0 & ! out , n2, s2, strig, ale_bl_stat & ! out , therm_tke_max, env_tke_max & ! out , alp_bl_det, alp_bl_fluct_m, alp_bl_fluct_tke & ! out , alp_bl_conv, alp_bl_stat & ! out ) USE indice_sol_mod USE lmdz_thermcell_main, ONLY: thermcell_tke_transport USE lmdz_alpale USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE INCLUDE "FCTTRE.h" !======================================================================= ! Auteurs: Catherine Rio ! Modifications : ! Nicolas Rochetin et Jean-Yves Grandpeix ! pour la fermeture stochastique. 2012 ! Frédéric Hourdin : ! netoyage informatique. 2022 !======================================================================= !----------------------------------------------------------------------- ! arguments: ! ---------- !------Entrees INTEGER, INTENT(IN) :: ngrid, nlay REAL, INTENT(IN) :: ptimestep REAL, INTENT(IN) :: pplay(ngrid, nlay), pplev(ngrid, nlay + 1) INTEGER, INTENT(IN), DIMENSION(ngrid) :: lmax, lalim REAL, INTENT(IN), DIMENSION(ngrid) :: zmax, zcong REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: zw2 REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: fraca REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: wth3 REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: rhobarz REAL, INTENT(IN), DIMENSION(ngrid) :: wmax_sec REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: entr0 REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: fm0, fm REAL, INTENT(IN), DIMENSION(ngrid) :: pcon REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: alim_star REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1, nbsrf) :: pbl_tke REAL, INTENT(IN), DIMENSION(ngrid, nbsrf) :: pctsrf REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: omega REAL, INTENT(IN), DIMENSION(ngrid) :: airephy !------Sorties REAL, INTENT(OUT), DIMENSION(ngrid) :: ale_bl, alp_bl REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: wght_th INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lalim_conv REAL, INTENT(OUT), DIMENSION(ngrid) :: zlcl, fraca0, w0, w_conv REAL, INTENT(OUT), DIMENSION(ngrid) :: therm_tke_max0, env_tke_max0, n2, s2, ale_bl_stat, strig REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: therm_tke_max, env_tke_max REAL, INTENT(OUT), DIMENSION(ngrid) :: alp_bl_det, alp_bl_fluct_m, alp_bl_fluct_tke REAL, INTENT(OUT), DIMENSION(ngrid) :: alp_bl_conv, alp_bl_stat !============================================================================================= !------Local !============================================================================================= REAL susqr2pi, reuler INTEGER ig, k, l INTEGER nsrf REAL rhobarz0(ngrid) ! Densité au LCL LOGICAL ok_lcl(ngrid) ! Existence du LCL des thermiques INTEGER klcl(ngrid) ! Niveau du LCL REAL interp(ngrid) ! Coef d'interpolation pour le LCL !--Triggering REAL, parameter :: su_cst = 4e4 ! Surface unite: celle d'un updraft élémentaire REAL, parameter :: hcoef = 1 ! Coefficient directeur pour le calcul de s2 REAL, parameter :: hmincoef = 0.3 ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 REAL, parameter :: eps1 = 0.3 ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd) REAL, DIMENSION(ngrid) :: hmin ! Ordonnée à l'origine pour le calcul de s2 REAL, DIMENSION(ngrid) :: zmax_moy ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl) REAL, parameter :: zmax_moy_coef = 0.33 REAL, DIMENSION(ngrid) :: depth ! Epaisseur moyenne du cumulus REAL, DIMENSION(ngrid) :: zcong_moy REAL, DIMENSION(ngrid) :: w_max ! Vitesse max statistique REAL, DIMENSION(ngrid) :: s_max(ngrid) !--Closure REAL, DIMENSION(ngrid, nlay) :: pbl_tke_max ! Profil de TKE moyenne REAL, DIMENSION(ngrid) :: pbl_tke_max0 ! TKE moyenne au LCL REAL, DIMENSION(ngrid, nlay) :: w_ls ! Vitesse verticale grande échelle (m/s) REAL, parameter :: coef_m = 1. ! On considère un rendement pour alp_bl_fluct_m REAL, parameter :: coef_tke = 1. ! On considère un rendement pour alp_bl_fluct_tke REAL :: zdp REAL, DIMENSION(ngrid) :: alp_int, dp_int REAL, DIMENSION(ngrid) :: fm_tot !------------------------------------------------------------ ! Initialize output arrays related to stochastic triggering !------------------------------------------------------------ DO ig = 1, ngrid zlcl(ig) = 0. fraca0(ig) = 0. w0(ig) = 0. w_conv(ig) = 0. therm_tke_max0(ig) = 0. env_tke_max0(ig) = 0. n2(ig) = 0. s2(ig) = 0. ale_bl_stat(ig) = 0. strig(ig) = 0. alp_bl_det(ig) = 0. alp_bl_fluct_m(ig) = 0. alp_bl_fluct_tke(ig) = 0. alp_bl_conv(ig) = 0. alp_bl_stat(ig) = 0. ENDDO DO l = 1, nlay DO ig = 1, ngrid therm_tke_max(ig, l) = 0. env_tke_max(ig, l) = 0. ENDDO ENDDO !------------Test sur le LCL des thermiques DO ig = 1, ngrid ok_lcl(ig) = .FALSE. IF ((pcon(ig) > pplay(ig, nlay - 1)) .AND. (pcon(ig) < pplay(ig, 1))) ok_lcl(ig) = .TRUE. enddo !------------Localisation des niveaux entourant le LCL et du coef d'interpolation DO l = 1, nlay - 1 DO ig = 1, ngrid IF (ok_lcl(ig)) THEN !ATTENTION,zw2 calcule en pplev ! if ((pplay(ig,l) .ge. pcon(ig)) .AND. (pplay(ig,l+1) .le. pcon(ig))) THEN ! klcl(ig)=l ! interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig))) ! endif IF ((pplev(ig, l) >= pcon(ig)) .AND. (pplev(ig, l + 1) <= pcon(ig))) THEN klcl(ig) = l interp(ig) = (pcon(ig) - pplev(ig, klcl(ig))) / (pplev(ig, klcl(ig) + 1) - pplev(ig, klcl(ig))) endif endif enddo enddo DO ig = 1, ngrid !CR:REHABILITATION ZMAX CONTINU IF (ok_lcl(ig)) THEN rhobarz0(ig) = rhobarz(ig, klcl(ig)) + (rhobarz(ig, klcl(ig) + 1) & - rhobarz(ig, klcl(ig))) * interp(ig) zlcl(ig) = (pplev(ig, 1) - pcon(ig)) / (rhobarz0(ig) * RG) zlcl(ig) = min(zlcl(ig), zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax else rhobarz0(ig) = 0. zlcl(ig) = zmax(ig) endif enddo !!jyg fin !------------Calcul des propriétés du thermique au LCL IF ((iflag_trig_bl>=1) .OR. (iflag_clos_bl>=1)) THEN !-----Initialisation de la TKE moyenne DO l = 1, nlay DO ig = 1, ngrid pbl_tke_max(ig, l) = 0. enddo enddo !-----Calcul de la TKE moyenne DO nsrf = 1, nbsrf DO l = 1, nlay DO ig = 1, ngrid pbl_tke_max(ig, l) = pctsrf(ig, nsrf) * pbl_tke(ig, l, nsrf) + pbl_tke_max(ig, l) enddo enddo enddo !-----Initialisations des TKE dans et hors des thermiques DO l = 1, nlay DO ig = 1, ngrid therm_tke_max(ig, l) = pbl_tke_max(ig, l) env_tke_max(ig, l) = pbl_tke_max(ig, l) enddo enddo !-----Calcul de la TKE transportée par les thermiques : therm_tke_max CALL thermcell_tke_transport(ngrid, nlay, ptimestep, fm0, entr0, & ! in rg, pplev, therm_tke_max) ! out ! PRINT *,' thermcell_tke_transport -> ' !!jyg !-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls DO l = 1, nlay DO ig = 1, ngrid pbl_tke_max(ig, l) = fraca(ig, l) * therm_tke_max(ig, l) + (1. - fraca(ig, l)) * env_tke_max(ig, l) ! Recalcul de TKE moyenne après transport de TKE_TH env_tke_max(ig, l) = (pbl_tke_max(ig, l) - fraca(ig, l) * therm_tke_max(ig, l)) / (1. - fraca(ig, l)) ! Recalcul de TKE dans l'environnement après transport de TKE_TH w_ls(ig, l) = -1. * omega(ig, l) / (RG * rhobarz(ig, l)) ! Vitesse verticale de grande échelle enddo enddo ! PRINT *,' apres w_ls = ' !!jyg DO ig = 1, ngrid IF (ok_lcl(ig)) THEN fraca0(ig) = fraca(ig, klcl(ig)) + (fraca(ig, klcl(ig) + 1) & - fraca(ig, klcl(ig))) * interp(ig) w0(ig) = zw2(ig, klcl(ig)) + (zw2(ig, klcl(ig) + 1) & - zw2(ig, klcl(ig))) * interp(ig) w_conv(ig) = w_ls(ig, klcl(ig)) + (w_ls(ig, klcl(ig) + 1) & - w_ls(ig, klcl(ig))) * interp(ig) therm_tke_max0(ig) = therm_tke_max(ig, klcl(ig)) & + (therm_tke_max(ig, klcl(ig) + 1) - therm_tke_max(ig, klcl(ig))) * interp(ig) env_tke_max0(ig) = env_tke_max(ig, klcl(ig)) + (env_tke_max(ig, klcl(ig) + 1) & - env_tke_max(ig, klcl(ig))) * interp(ig) pbl_tke_max0(ig) = pbl_tke_max(ig, klcl(ig)) + (pbl_tke_max(ig, klcl(ig) + 1) & - pbl_tke_max(ig, klcl(ig))) * interp(ig) IF (therm_tke_max0(ig)>=20.) therm_tke_max0(ig) = 20. IF (env_tke_max0(ig)>=20.) env_tke_max0(ig) = 20. IF (pbl_tke_max0(ig)>=20.) pbl_tke_max0(ig) = 20. else fraca0(ig) = 0. w0(ig) = 0. !!jyg le 27/04/2012 !! zlcl(ig)=0. !! endif enddo ENDIF ! IF ( (iflag_trig_bl.ge.1) .OR. (iflag_clos_bl.ge.1) ) ! PRINT *,'ENDIF ( (iflag_trig_bl.ge.1) .OR. (iflag_clos_bl.ge.1) ) ' !!jyg !------------Triggering------------------ IF (iflag_trig_bl>=1) THEN !-----Initialisations depth(:) = 0. n2(:) = 0. s2(:) = 100. ! some low value, arbitrary s_max(:) = 0. !-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max) DO ig = 1, ngrid zmax_moy(ig) = zlcl(ig) + zmax_moy_coef * (zmax(ig) - zlcl(ig)) depth(ig) = zmax_moy(ig) - zlcl(ig) hmin(ig) = hmincoef * zlcl(ig) IF (depth(ig)>=10.) THEN s2(ig) = (hcoef * depth(ig) + hmin(ig))**2 n2(ig) = (1. - eps1) * fraca0(ig) * airephy(ig) / s2(ig) !! !!jyg le 27/04/2012 !! s_max(ig)=s2(ig)*log(n2(ig)) !! if (n2(ig) .lt. 1) s_max(ig)=0. s_max(ig) = s2(ig) * log(max(n2(ig), 1.)) !!fin jyg else n2(ig) = 0. s_max(ig) = 0. endif enddo ! PRINT *,'avant Calcul de Wmax ' !!jyg !CR: calcul de strig IF (iflag_strig==0) THEN strig(:) = s_trig ELSE IF (iflag_strig==1) THEN DO ig = 1, ngrid ! zcong_moy(ig)=zlcl(ig)+zmax_moy_coef*(zcong(ig)-zlcl(ig)) ! strig(ig)=(hcoef*(zcong_moy(ig)-zlcl(ig))+hmin(ig))**2 strig(ig) = (zcong(ig) - zlcl(ig))**2 enddo ELSE IF (iflag_strig==2) THEN DO ig = 1, ngrid IF (h_trig>zlcl(ig)) THEN strig(ig) = (h_trig - zlcl(ig))**2 else strig(ig) = s_trig endif enddo endif susqr2pi = su_cst * sqrt(2. * Rpi) reuler = exp(1.) DO ig = 1, ngrid IF ((depth(ig)>=10.) .AND. (s_max(ig)>susqr2pi * reuler)) THEN w_max(ig) = w0(ig) * (1. + sqrt(2. * log(s_max(ig) / susqr2pi) - log(2. * log(s_max(ig) / susqr2pi)))) ale_bl_stat(ig) = 0.5 * w_max(ig)**2 else w_max(ig) = 0. ale_bl_stat(ig) = 0. endif enddo ENDIF ! iflag_trig_bl ! PRINT *,'ENDIF iflag_trig_bl' !!jyg !------------Closure------------------ IF (iflag_clos_bl>=2) THEN !-----Calcul de ALP_BL_STAT DO ig = 1, ngrid alp_bl_det(ig) = 0.5 * coef_m * rhobarz0(ig) * (w0(ig)**3) * fraca0(ig) * (1. - 2. * fraca0(ig)) / ((1. - fraca0(ig))**2) alp_bl_fluct_m(ig) = 1.5 * rhobarz0(ig) * fraca0(ig) * (w_conv(ig) + coef_m * w0(ig)) * & (w0(ig)**2) alp_bl_fluct_tke(ig) = 3. * coef_m * rhobarz0(ig) * w0(ig) * fraca0(ig) * (therm_tke_max0(ig) - env_tke_max0(ig)) & + 3. * rhobarz0(ig) * w_conv(ig) * pbl_tke_max0(ig) IF (iflag_clos_bl>=2) THEN alp_bl_conv(ig) = 1.5 * coef_m * rhobarz0(ig) * fraca0(ig) * (fraca0(ig) / (1. - fraca0(ig))) * w_conv(ig) * & (w0(ig)**2) else alp_bl_conv(ig) = 0. endif alp_bl_stat(ig) = alp_bl_det(ig) + alp_bl_fluct_m(ig) + alp_bl_fluct_tke(ig) + alp_bl_conv(ig) enddo !-----Sécurité ALP infinie DO ig = 1, ngrid IF (fraca0(ig)>0.98) alp_bl_stat(ig) = 2. enddo ENDIF ! (iflag_clos_bl.ge.2) !!! fin nrlmd le 10/04/2012 ! PRINT*,'avant calcul ale et alp' !calcul de ALE et ALP pour la convection alp_bl(:) = 0. ale_bl(:) = 0. ! PRINT*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)' DO l = 1, nlay DO ig = 1, ngrid alp_bl(ig) = max(alp_bl(ig), 0.5 * rhobarz(ig, l) * wth3(ig, l)) ale_bl(ig) = max(ale_bl(ig), 0.5 * zw2(ig, l)**2) ! PRINT*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig) enddo enddo ! ale sec (max de wmax/2 sous la zone d'inhibition) dans ! le cas iflag_trig_bl=3 IF (iflag_trig_bl==3) ale_bl(:) = 0.5 * wmax_sec(:)**2 !test:calcul de la ponderation des couches pour KE !initialisations fm_tot(:) = 0. wght_th(:, :) = 1. lalim_conv(:) = lalim(:) DO k = 1, nlay DO ig = 1, ngrid IF (k<=lalim_conv(ig)) fm_tot(ig) = fm_tot(ig) + fm(ig, k) enddo enddo ! assez bizarre car, si on est dans la couche d'alim et que alim_star et ! plus petit que 1.e-10, on prend wght_th=1. DO k = 1, nlay DO ig = 1, ngrid IF (k<=lalim_conv(ig).AND.alim_star(ig, k)>1.e-10) THEN wght_th(ig, k) = alim_star(ig, k) 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)<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 !------------------------------------------------------------------------ ! Modif CR/FH 20110310 : alp integree sur la verticale. ! Integrale verticale de ALP. ! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des ! couches !------------------------------------------------------------------------ alp_int(:) = 0. dp_int(:) = 0. DO l = 2, nlay DO ig = 1, ngrid IF(l<=lmax(ig)) THEN zdp = pplay(ig, l - 1) - pplay(ig, l) alp_int(ig) = alp_int(ig) + 0.5 * rhobarz(ig, l) * wth3(ig, l) * zdp dp_int(ig) = dp_int(ig) + zdp endif enddo enddo IF (iflag_coupl>=3 .AND. iflag_coupl<=5) THEN DO ig = 1, ngrid !valeur integree de alp_bl * 0.5: IF (dp_int(ig)>0.) THEN alp_bl(ig) = alp_int(ig) / dp_int(ig) endif enddo! endif ! Facteur multiplicatif sur alp_bl alp_bl(:) = alp_bl_k * alp_bl(:) !------------------------------------------------------------------------ RETURN END END MODULE lmdz_thermcell_alp