Context Navigation

alpale_th.F90 @ 4826

Last change on this file since 4826 was 4826, checked in by jyg, 3 months ago
New flag in alpale_th.F90: multiply_proba_notrig; default=.FALSE.
Property svn:keywords set to `Id`
File size: 12.6 KB

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

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

Download in other formats:

Original Format