Context Navigation

alpale_th.F90 @ 3605

Last change on this file since 3605 was 3531, checked in by Laurent Fairhead, 5 years ago
Replaced STOP statements by a call to abort_physic in phylmd as per ticket #86 Still some work to be done in phylmd subdirectories
Property svn:keywords set to `Id`
File size: 12.2 KB

Rev	Line
[3531]	1	!
	2	! $Id: alpale_th.F90 3531 2019-06-06 15:08:45Z lguez $
	3	!
[2558]	4	SUBROUTINE alpale_th ( dtime, lmax_th, t_seri, cell_area, &
[2513]	5	cin, s2, n2, &
	6	ale_bl_trig, ale_bl_stat, ale_bl, &
[2556]	7	alp_bl, alp_bl_stat, &
[3209]	8	proba_notrig, random_notrig, birth_rate)
[2513]	9
	10	! **************************************************************
	11	! *
	12	! ALPALE_TH *
	13	! *
	14	! *
	15	! written by : Jean-Yves Grandpeix, 11/05/2016 *
	16	! modified by : *
	17	! **************************************************************
	18
	19	USE dimphy
	20	USE ioipsl_getin_p_mod, ONLY : getin_p
	21	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
	22	!
	23	IMPLICIT NONE
	24
	25	!================================================================
	26	! Auteur(s) : Jean-Yves Grandpeix, 11/05/2016
	27	! Objet : Contribution of the thermal scheme to Ale and Alp
	28	!================================================================
	29
	30	! Input arguments
	31	!----------------
	32	REAL, INTENT(IN) :: dtime
[2558]	33	REAL, DIMENSION(klon), INTENT(IN) :: cell_area
[2513]	34	INTEGER, DIMENSION(klon), INTENT(IN) :: lmax_th
	35	REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri
	36	REAL, DIMENSION(klon), INTENT(IN) :: ale_bl_stat
	37	REAL, DIMENSION(klon), INTENT(IN) :: cin
	38	REAL, DIMENSION(klon), INTENT(IN) :: s2, n2
	39
	40	REAL, DIMENSION(klon), INTENT(INOUT) :: ale_bl_trig, ale_bl
	41	REAL, DIMENSION(klon), INTENT(INOUT) :: alp_bl
	42	REAL, DIMENSION(klon), INTENT(INOUT) :: alp_bl_stat
[2730]	43	REAL, DIMENSION(klon), INTENT(INOUT) :: proba_notrig
[2513]	44
[2730]	45	REAL, DIMENSION(klon), INTENT(OUT) :: random_notrig
[2556]	46
[3208]	47	REAL, DIMENSION(klon), INTENT(OUT) :: birth_rate
	48
[2513]	49	include "thermcell.h"
	50
	51	! Local variables
	52	!----------------
	53	INTEGER :: i
	54	LOGICAL, SAVE :: first = .TRUE.
	55	REAL, SAVE :: random_notrig_max=1.
[2558]	56	REAL, SAVE :: cv_feed_area
	57	REAL :: birth_number
	58	REAL, DIMENSION(klon) :: ale_bl_ref
[2556]	59	REAL, DIMENSION(klon) :: tau_trig
[2513]	60	!
	61	!$OMP THREADPRIVATE(random_notrig_max)
[2558]	62	!$OMP THREADPRIVATE(cv_feed_area)
[2513]	63	!$OMP THREADPRIVATE(first)
	64	!
[2558]	65	REAL umexp ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
	66	REAL x
[3531]	67	CHARACTER (LEN=20) :: modname='alpale_th'
	68	CHARACTER (LEN=80) :: abort_message
	69
[2558]	70	umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
[3073]	71	(1.-max(sign(1.,x-1.e-3),0.))(1.-0.5x(1.-x/3.(1.-0.25*x))) !!! correct formula (jyg)
	72	!!! (1.-max(sign(1.,x-1.e-3),0.))(-0.5x(1.-x/3.(1.-0.25*x))) !!! bug introduced by mistake (jyg)
	73	!!! (1.-max(sign(1.,x-1.e-3),0.))(1.-0.5x(1.-x/3.(1.-0.25*x))) !!! initial correct formula (jyg)
[2558]	74	!
	75	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	76	! JYG, 20160513 : Introduction of the Effective Lifting Power (ELP), which
	77	! takes into account the area (cv_feed_area) covered by thermals contributing
	78	! to each cumulonimbus.
	79	! The use of ELP prevents singularities when the trigger probability tends to
	80	! zero. It is activated by iflag_clos_bl = 3.
	81	! The ELP values are stored in the ALP_bl variable.
	82	!
	83	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	84	!
[2730]	85	!!
	86	!! The following 3 lines should be commented if one wants to activate the
	87	!! multiplication of no-trigger probabilities between calls to the convection
	88	!! scheme.
	89	!!
	90	do i=1,klon
	91	proba_notrig(i)=1.
	92	enddo
	93	!!
	94	!!
[2558]	95	!---------------------------------------
	96	IF (iflag_clos_bl .LT. 3) THEN
	97	!---------------------------------------
	98	!
	99	! Original code (Nicolas Rochetin)
	100	! --------------------------------
	101
[2513]	102	IF (first) THEN
[2558]	103	random_notrig_max=1.
[2513]	104	CALL getin_p('random_notrig_max',random_notrig_max)
[2515]	105	first=.FALSE.
[2513]	106	ENDIF
	107	!cc nrlmd le 10/04/2012
	108	!-----------Stochastic triggering-----------
	109	if (iflag_trig_bl.ge.1) then
	110	!
	111	IF (prt_level .GE. 10) THEN
[3531]	112	WRITE(lunout,*)'cin, ale_bl_stat, alp_bl, alp_bl_stat ', &
[3208]	113	cin, ale_bl_stat, alp_bl, alp_bl_stat
[2513]	114	ENDIF
	115
	116
	117	!----Initialisations
	118	do i=1,klon
[2730]	119	!!jyg proba_notrig(i)=1.
[2513]	120	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
	121	if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
	122	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
	123	tau_trig(i)=tau_trig_shallow
	124	else
	125	tau_trig(i)=tau_trig_deep
	126	endif
	127	enddo
	128	!
	129	IF (prt_level .GE. 10) THEN
[3531]	130	WRITE(lunout,*)'random_notrig, tau_trig ', &
[2513]	131	random_notrig, tau_trig
[3531]	132	WRITE(lunout,*)'s_trig,s2,n2 ', &
[2513]	133	s_trig,s2,n2
	134	ENDIF
	135
	136	!Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
	137	IF (iflag_trig_bl.eq.1) then
	138
	139	!----Tirage al\'eatoire et calcul de ale_bl_trig
	140	do i=1,klon
	141	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
[2730]	142	proba_notrig(i)=proba_notrig(i)* &
	143	(1.-exp(-s_trig/s2(i)))*(n2(i)dtime/tau_trig(i))
[2513]	144	! print *, 'proba_notrig(i) ',proba_notrig(i)
	145	if (random_notrig(i) .ge. proba_notrig(i)) then
	146	ale_bl_trig(i)=ale_bl_stat(i)
	147	else
	148	ale_bl_trig(i)=0.
	149	endif
[3208]	150	birth_rate(i) = n2(i)exp(-s_trig/s2(i))/(tau_trig(i)cell_area(i))
	151	!!! birth_rate(i) = max(birth_rate(i),1.e-18)
[2513]	152	else
[2730]	153	!!jyg proba_notrig(i)=1.
[3208]	154	birth_rate(i) = 0.
[2513]	155	random_notrig(i)=0.
	156	ale_bl_trig(i)=0.
	157	endif
	158	enddo
	159
	160	ELSE IF (iflag_trig_bl.ge.2) then
	161
	162	do i=1,klon
	163	if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) ) then
[2730]	164	proba_notrig(i)=proba_notrig(i)* &
	165	(1.-exp(-s_trig/s2(i)))*(n2(i)dtime/tau_trig(i))
[2513]	166	! print *, 'proba_notrig(i) ',proba_notrig(i)
	167	if (random_notrig(i) .ge. proba_notrig(i)) then
	168	ale_bl_trig(i)=Ale_bl(i)
	169	else
	170	ale_bl_trig(i)=0.
	171	endif
[3208]	172	birth_rate(i) = n2(i)exp(-s_trig/s2(i))/(tau_trig(i)cell_area(i))
	173	!!! birth_rate(i) = max(birth_rate(i),1.e-18)
[2513]	174	else
[2730]	175	!!jyg proba_notrig(i)=1.
[3208]	176	birth_rate(i) = 0.
[2513]	177	random_notrig(i)=0.
	178	ale_bl_trig(i)=0.
	179	endif
	180	enddo
	181
	182	ENDIF
	183
	184	!
	185	IF (prt_level .GE. 10) THEN
[3531]	186	WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
[2513]	187	proba_notrig, ale_bl_trig
	188	ENDIF
	189
	190	endif !(iflag_trig_bl)
	191
	192	!-----------Statistical closure-----------
	193	if (iflag_clos_bl.eq.1) then
	194
	195	do i=1,klon
	196	!CR: alp probabiliste
	197	if (ale_bl_trig(i).gt.0.) then
	198	alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
	199	endif
	200	enddo
	201
	202	else if (iflag_clos_bl.eq.2) then
	203
	204	!CR: alp calculee dans thermcell_main
	205	do i=1,klon
	206	alp_bl(i)=alp_bl_stat(i)
	207	enddo
	208
	209	else
	210
	211	alp_bl_stat(:)=0.
	212
	213	endif !(iflag_clos_bl)
	214
[2558]	215	!
	216	!---------------------------------------
	217	ELSEIF (iflag_clos_bl .EQ. 3) THEN ! (iflag_clos_bl .LT. 3)
	218	!---------------------------------------
	219	!
	220	! New code with Effective Lifting Power
	221	! -------------------------------------
	222	IF (first) THEN
	223	cv_feed_area = 1.e10 ! m2
	224	CALL getin_p('cv_feed_area', cv_feed_area)
	225	first=.FALSE.
	226	ENDIF
	227
	228	!-----------Stochastic triggering-----------
	229	if (iflag_trig_bl.ge.1) then
	230	!
	231	IF (prt_level .GE. 10) THEN
[3531]	232	WRITE(lunout,*)'cin, ale_bl_stat, alp_bl_stat ', &
[2558]	233	cin, ale_bl_stat, alp_bl_stat
	234	ENDIF
	235
	236	! Use ale_bl_stat (Rochetin's code) or ale_bl (old code) according to
	237	! iflag_trig_bl value.
	238	IF (iflag_trig_bl.eq.1) then ! use ale_bl_stat (Rochetin computation)
	239	do i=1,klon
	240	ale_bl_ref(i)=ale_bl_stat(i)
	241	enddo
	242	ELSE IF (iflag_trig_bl.ge.2) then ! use ale_bl (old computation)
	243	do i=1,klon
	244	ale_bl_ref(i)=Ale_bl(i)
	245	enddo
	246	ENDIF ! (iflag_trig_bl.eq.1)
	247
	248
	249	!----Initializations and random number generation
	250	do i=1,klon
[2730]	251	!!jyg proba_notrig(i)=1.
[2558]	252	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
	253	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
	254	tau_trig(i)=tau_trig_shallow
	255	else
	256	tau_trig(i)=tau_trig_deep
	257	endif
	258	enddo
	259	!
	260	IF (prt_level .GE. 10) THEN
[3531]	261	WRITE(lunout,*)'random_notrig, tau_trig ', &
[2558]	262	random_notrig, tau_trig
[3531]	263	WRITE(lunout,*)'s_trig,s2,n2 ', &
[2558]	264	s_trig,s2,n2
	265	ENDIF
	266
	267	!----alp_bl computation
	268	do i=1,klon
	269	if ( (ale_bl_ref(i) .gt. abs(cin(i))+1.e-10) ) then
	270	birth_number = n2(i)*exp(-s_trig/s2(i))
	271	birth_rate(i) = birth_number/(tau_trig(i)*cell_area(i))
[3208]	272	!!! birth_rate(i) = max(birth_rate(i),1.e-18)
[2730]	273	proba_notrig(i)=proba_notrig(i)exp(-birth_numberdtime/tau_trig(i))
[2558]	274	Alp_bl(i) = Alp_bl(i)* &
	275	umexp(-birth_number*cv_feed_area/cell_area(i))/ &
	276	umexp(-birth_numberdtime/tau_trig(i)) &
	277	tau_trig(i)cv_feed_area/(dtimecell_area(i))
	278	else
[2730]	279	!!jyg proba_notrig(i)=1.
[3208]	280	birth_rate(i)=0.
[2558]	281	random_notrig(i)=0.
	282	alp_bl(i)=0.
	283	endif
	284	enddo
	285
	286	!----ale_bl_trig computation
	287	do i=1,klon
	288	if (random_notrig(i) .ge. proba_notrig(i)) then
	289	ale_bl_trig(i)=ale_bl_ref(i)
	290	else
	291	ale_bl_trig(i)=0.
	292	endif
	293	enddo
	294
	295	!
	296	IF (prt_level .GE. 10) THEN
[3531]	297	WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
[2558]	298	proba_notrig, ale_bl_trig
	299	ENDIF
	300
	301	endif !(iflag_trig_bl .ge. 1)
	302
	303	!---------------------------------------
	304	ENDIF ! (iflag_clos_bl .LT. 3)
	305	!---------------------------------------
	306
[2513]	307	IF (prt_level .GE. 10) THEN
[3531]	308	WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', &
[3209]	309	ale_bl_trig(1), alp_bl_stat(1), birth_rate(1)
[2513]	310	ENDIF
	311
	312	!cc fin nrlmd le 10/04/2012
[2554]	313	!
[2513]	314	!IM/FH: 2011/02/23
	315	! Couplage Thermiques/Emanuel seulement si T<0
	316	if (iflag_coupl==2) then
	317	IF (prt_level .GE. 10) THEN
[3531]	318	WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0'
[2513]	319	ENDIF
	320	do i=1,klon
	321	if (t_seri(i,lmax_th(i))>273.) then
	322	Ale_bl(i)=0.
	323	endif
	324	enddo
[3531]	325	! print *,'In order to run with iflag_coupl=2, you have to comment out the following stop'
	326	! STOP
	327	abort_message='In order to run with iflag_coupl=2, you have to comment out the following abort'
	328	CALL abort_physic(modname,abort_message,1)
[2513]	329	endif
	330	RETURN
	331	END
	332

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format