! $Id: alpale_th.F90 5160 2024-08-03 12:56:58Z evignon $ SUBROUTINE alpale_th( dtime, lmax_th, t_seri, cell_area, & cin, s2, n2, strig, & ale_bl_trig, ale_bl_stat, ale_bl, & alp_bl, alp_bl_stat, & proba_notrig, random_notrig, birth_rate) ! ************************************************************** ! * ! ALPALE_TH * ! * ! * ! written by : Jean-Yves Grandpeix, 11/05/2016 * ! modified by : * ! ************************************************************** USE dimphy USE lmdz_ioipsl_getin_p, ONLY: getin_p USE lmdz_print_control, ONLY: mydebug=>debug , lunout, prt_level USE lmdz_abort_physic, ONLY: abort_physic USE lmdz_alpale USE lmdz_cv, ONLY: cv_feed IMPLICIT NONE !================================================================ ! Auteur(s) : Jean-Yves Grandpeix, 11/05/2016 ! Objet : Contribution of the thermal scheme to Ale and Alp !================================================================ ! Input arguments !---------------- REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(IN) :: cell_area INTEGER, DIMENSION(klon), INTENT(IN) :: lmax_th REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri REAL, DIMENSION(klon), INTENT(IN) :: ale_bl_stat REAL, DIMENSION(klon), INTENT(IN) :: cin REAL, DIMENSION(klon), INTENT(IN) :: s2, n2, strig REAL, DIMENSION(klon), INTENT(INOUT) :: ale_bl_trig, ale_bl REAL, DIMENSION(klon), INTENT(INOUT) :: alp_bl REAL, DIMENSION(klon), INTENT(INOUT) :: alp_bl_stat REAL, DIMENSION(klon), INTENT(INOUT) :: proba_notrig REAL, DIMENSION(klon), INTENT(OUT) :: random_notrig REAL, DIMENSION(klon), INTENT(OUT) :: birth_rate ! Local variables !---------------- INTEGER :: i LOGICAL, SAVE :: first = .TRUE. LOGICAL, SAVE :: multiply_proba_notrig = .FALSE. REAL, SAVE :: random_notrig_max=1. REAL, SAVE :: cv_feed_area REAL :: birth_number REAL, DIMENSION(klon) :: ale_bl_ref REAL, DIMENSION(klon) :: tau_trig !$OMP THREADPRIVATE(multiply_proba_notrig) !$OMP THREADPRIVATE(random_notrig_max) !$OMP THREADPRIVATE(cv_feed_area) !$OMP THREADPRIVATE(first) REAL umexp ! expression of (1.-exp(-x))/x valid for all x, especially when x->0 REAL x CHARACTER (LEN=20) :: modname='alpale_th' CHARACTER (LEN=80) :: abort_message umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + & (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! correct formula (jyg) !!! (1.-max(sign(1.,x-1.e-3),0.))*(-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! bug introduced by mistake (jyg) !!! (1.-max(sign(1.,x-1.e-3),0.))*(1.-0.5*x*(1.-x/3.*(1.-0.25*x))) !!! initial correct formula (jyg) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which ! takes into account the area (cv_feed_area) covered by thermals contributing ! to each cumulonimbus. ! The use of ELP prevents singularities when the trigger probability tends to ! zero. It is activated by iflag_clos_bl = 3. ! The ELP values are stored in the ALP_bl variable. !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF (first) THEN CALL getin_p('multiply_proba_notrig',multiply_proba_notrig) IF (iflag_clos_bl < 3) THEN random_notrig_max=1. CALL getin_p('random_notrig_max',random_notrig_max) ELSEIF (iflag_clos_bl == 3) THEN ! (iflag_clos_bl .LT. 3) cv_feed_area = 1.e10 ! m2 CALL getin_p('cv_feed_area', cv_feed_area) ENDIF !! (iflag_clos_bl .LT. 3) first=.FALSE. ENDIF !! !! Control of the multiplication of no-trigger probabilities between calls !! to the convection scheme. If multiply_proba_notrig is .FALSE., THEN !! proba_notrig is set to 1 at each CALL to alpale_th, so that only the last CALL !! plays a role in the triggering of convection. If it is .TRUE., then propa_notrig !! is saved between calls to convection and is reset to 1 only after calling the !! convection scheme. !! For instance, if the probability of no_trigger is 0.9 at each call, and if !! there are 3 calls to alpale_th between calls to the convection scheme, then the !! probability of triggering convection will be 0.1 (= 1.-0.9) if !! multiply_proba_notrig is .FALSE. and 0.271 (= 1.-0.9^3) if multiply_proba_notrig !! is .TRUE. !! IF (.NOT.multiply_proba_notrig) THEN DO i=1,klon proba_notrig(i)=1. ENDDO ENDIF !! (.NOT.multiply_proba_notrig) !! !! !--------------------------------------- IF (iflag_clos_bl < 3) THEN !--------------------------------------- ! Original code (Nicolas Rochetin) ! -------------------------------- !cc nrlmd le 10/04/2012 !-----------Stochastic triggering----------- IF (iflag_trig_bl>=1) THEN IF (prt_level >= 10) THEN WRITE(lunout,*)'cin, ale_bl_stat, alp_bl, alp_bl_stat ', & cin, ale_bl_stat, alp_bl, alp_bl_stat ENDIF !----Initialisations DO i=1,klon !!jyg proba_notrig(i)=1. random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i)) IF ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0. IF ( ale_bl_trig(i) < abs(cin(i))+1.e-10 ) THEN tau_trig(i)=tau_trig_shallow else tau_trig(i)=tau_trig_deep endif enddo IF (prt_level >= 10) THEN WRITE(lunout,*)'random_notrig, tau_trig ', & random_notrig, tau_trig WRITE(lunout,*)'s_trig,s2,n2 ', & s_trig,s2,n2 ENDIF !Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2) IF (iflag_trig_bl==1) THEN !----Tirage al\'eatoire et calcul de ale_bl_trig DO i=1,klon IF ( (ale_bl_stat(i) > abs(cin(i))+1.e-10) ) THEN proba_notrig(i)=proba_notrig(i)* & (1.-exp(-strig(i)/s2(i)))**(n2(i)*dtime/tau_trig(i)) ! PRINT *, 'proba_notrig(i) ',proba_notrig(i) IF (random_notrig(i) >= proba_notrig(i)) THEN ale_bl_trig(i)=ale_bl_stat(i) else ale_bl_trig(i)=0. endif birth_rate(i) = n2(i)*exp(-strig(i)/s2(i))/(tau_trig(i)*cell_area(i)) !!! birth_rate(i) = max(birth_rate(i),1.e-18) else !!jyg proba_notrig(i)=1. birth_rate(i) = 0. random_notrig(i)=0. ale_bl_trig(i)=0. endif enddo ELSE IF (iflag_trig_bl>=2) THEN DO i=1,klon IF ( (Ale_bl(i) > abs(cin(i))+1.e-10) ) THEN proba_notrig(i)=proba_notrig(i)* & (1.-exp(-strig(i)/s2(i)))**(n2(i)*dtime/tau_trig(i)) ! PRINT *, 'proba_notrig(i) ',proba_notrig(i) IF (random_notrig(i) >= proba_notrig(i)) THEN ale_bl_trig(i)=Ale_bl(i) else ale_bl_trig(i)=0. endif birth_rate(i) = n2(i)*exp(-strig(i)/s2(i))/(tau_trig(i)*cell_area(i)) !!! birth_rate(i) = max(birth_rate(i),1.e-18) else !!jyg proba_notrig(i)=1. birth_rate(i) = 0. random_notrig(i)=0. ale_bl_trig(i)=0. endif enddo ENDIF IF (prt_level >= 10) THEN WRITE(lunout,*)'proba_notrig, ale_bl_trig ', & proba_notrig, ale_bl_trig ENDIF endif !(iflag_trig_bl) !-----------Statistical closure----------- IF (iflag_clos_bl==1) THEN DO i=1,klon !CR: alp probabiliste IF (ale_bl_trig(i)>0.) THEN alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999)) endif enddo ELSE IF (iflag_clos_bl==2) THEN !CR: alp calculee dans thermcell_main DO i=1,klon alp_bl(i)=alp_bl_stat(i) enddo else alp_bl_stat(:)=0. endif !(iflag_clos_bl) !--------------------------------------- ELSEIF (iflag_clos_bl == 3) THEN ! (iflag_clos_bl .LT. 3) !--------------------------------------- ! New code with Effective Lifting Power ! ------------------------------------- !-----------Stochastic triggering----------- IF (iflag_trig_bl>=1) THEN IF (prt_level >= 10) THEN WRITE(lunout,*)'cin, ale_bl_stat, alp_bl_stat ', & cin, ale_bl_stat, alp_bl_stat ENDIF ! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to ! iflag_trig_bl value. IF (iflag_trig_bl==1) then ! use ale_bl_stat (Rochetin computation) DO i=1,klon ale_bl_ref(i)=ale_bl_stat(i) enddo ELSE IF (iflag_trig_bl>=2) then ! use ale_bl (old computation) DO i=1,klon ale_bl_ref(i)=Ale_bl(i) enddo ENDIF ! (iflag_trig_bl.EQ.1) !----Initializations and random number generation DO i=1,klon !!jyg proba_notrig(i)=1. random_notrig(i)=1e6*ale_bl_stat(i)-int(1e6*ale_bl_stat(i)) IF ( ale_bl_trig(i) < abs(cin(i))+1.e-10 ) THEN tau_trig(i)=tau_trig_shallow else tau_trig(i)=tau_trig_deep endif enddo IF (prt_level >= 10) THEN WRITE(lunout,*)'random_notrig, tau_trig ', & random_notrig, tau_trig WRITE(lunout,*)'s_trig,s2,n2 ', & s_trig,s2,n2 ENDIF !----alp_bl computation DO i=1,klon IF ( (ale_bl_ref(i) > abs(cin(i))+1.e-10) ) THEN birth_number = n2(i)*exp(-strig(i)/s2(i)) birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i)) !!! birth_rate(i) = max(birth_rate(i),1.e-18) proba_notrig(i)=proba_notrig(i)*exp(-birth_number*dtime/tau_trig(i)) Alp_bl(i) = Alp_bl(i)* & umexp(-birth_number*cv_feed_area/cell_area(i))/ & umexp(-birth_number*dtime/tau_trig(i))* & tau_trig(i)*cv_feed_area/(dtime*cell_area(i)) else !!jyg proba_notrig(i)=1. birth_rate(i)=0. random_notrig(i)=0. alp_bl(i)=0. endif enddo !----ale_bl_trig computation DO i=1,klon IF (random_notrig(i) >= proba_notrig(i)) THEN ale_bl_trig(i)=ale_bl_ref(i) else ale_bl_trig(i)=0. endif enddo IF (prt_level >= 10) THEN WRITE(lunout,*)'proba_notrig, ale_bl_trig ', & proba_notrig, ale_bl_trig ENDIF endif !(iflag_trig_bl .ge. 1) !--------------------------------------- ENDIF ! (iflag_clos_bl .LT. 3) !--------------------------------------- IF (prt_level >= 10) THEN WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', & ale_bl_trig(1), alp_bl_stat(1), birth_rate(1) ENDIF !cc fin nrlmd le 10/04/2012 !IM/FH: 2011/02/23 ! Couplage Thermiques/Emanuel seulement si T<0 IF (iflag_coupl==2) THEN IF (prt_level >= 10) THEN WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0' ENDIF DO i=1,klon IF (t_seri(i,lmax_th(i))>273.) THEN Ale_bl(i)=0. endif enddo ! PRINT *,'In order to run with iflag_coupl=2, you have to comment out the following stop ' ! STOP abort_message='In order to run with iflag_coupl=2, you have to comment out the following abort' CALL abort_physic(modname,abort_message,1) endif RETURN END