Context Navigation

alpale_th.F90 @ 5029

Last change on this file since 5029 was 3073, checked in by jyg, 8 years ago
back to original (correct) formula in alpale_th
File size: 11.3 KB

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

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

Download in other formats:

Original Format