Context Navigation

alpale_th.F90 @ 5321

Last change on this file since 5321 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

Line
1	!
2	! $Id: alpale_th.F90 3531 2019-06-06 15:08:45Z fairhead $
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 "thermcell.h"
50
51	! Local variables
52	!----------------
53	INTEGER :: i
54	LOGICAL, SAVE :: first = .TRUE.
55	REAL, SAVE :: random_notrig_max=1.
56	REAL, SAVE :: cv_feed_area
57	REAL :: birth_number
58	REAL, DIMENSION(klon) :: ale_bl_ref
59	REAL, DIMENSION(klon) :: tau_trig
60	!
61	!$OMP THREADPRIVATE(random_notrig_max)
62	!$OMP THREADPRIVATE(cv_feed_area)
63	!$OMP THREADPRIVATE(first)
64	!
65	REAL umexp ! expression of (1.-exp(-x))/x valid for all x, especially when x->0
66	REAL x
67	CHARACTER (LEN=20) :: modname='alpale_th'
68	CHARACTER (LEN=80) :: abort_message
69
70	umexp(x) = max(sign(1.,x-1.e-3),0.)*(1.-exp(-x))/max(x,1.e-3) + &
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)
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	!
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	!!
95	!---------------------------------------
96	IF (iflag_clos_bl .LT. 3) THEN
97	!---------------------------------------
98	!
99	! Original code (Nicolas Rochetin)
100	! --------------------------------
101
102	IF (first) THEN
103	random_notrig_max=1.
104	CALL getin_p('random_notrig_max',random_notrig_max)
105	first=.FALSE.
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
112	WRITE(lunout,*)'cin, ale_bl_stat, alp_bl, alp_bl_stat ', &
113	cin, ale_bl_stat, alp_bl, alp_bl_stat
114	ENDIF
115
116
117	!----Initialisations
118	do i=1,klon
119	!!jyg proba_notrig(i)=1.
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
130	WRITE(lunout,*)'random_notrig, tau_trig ', &
131	random_notrig, tau_trig
132	WRITE(lunout,*)'s_trig,s2,n2 ', &
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
142	proba_notrig(i)=proba_notrig(i)* &
143	(1.-exp(-s_trig/s2(i)))*(n2(i)dtime/tau_trig(i))
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
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)
152	else
153	!!jyg proba_notrig(i)=1.
154	birth_rate(i) = 0.
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
164	proba_notrig(i)=proba_notrig(i)* &
165	(1.-exp(-s_trig/s2(i)))*(n2(i)dtime/tau_trig(i))
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
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)
174	else
175	!!jyg proba_notrig(i)=1.
176	birth_rate(i) = 0.
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
186	WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
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
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
232	WRITE(lunout,*)'cin, ale_bl_stat, alp_bl_stat ', &
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
251	!!jyg proba_notrig(i)=1.
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
261	WRITE(lunout,*)'random_notrig, tau_trig ', &
262	random_notrig, tau_trig
263	WRITE(lunout,*)'s_trig,s2,n2 ', &
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))
272	!!! birth_rate(i) = max(birth_rate(i),1.e-18)
273	proba_notrig(i)=proba_notrig(i)exp(-birth_numberdtime/tau_trig(i))
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
279	!!jyg proba_notrig(i)=1.
280	birth_rate(i)=0.
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
297	WRITE(lunout,*)'proba_notrig, ale_bl_trig ', &
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
307	IF (prt_level .GE. 10) THEN
308	WRITE(lunout,*)'alpale_th: ale_bl_trig, alp_bl_stat, birth_rate ', &
309	ale_bl_trig(1), alp_bl_stat(1), birth_rate(1)
310	ENDIF
311
312	!cc fin nrlmd le 10/04/2012
313	!
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
318	WRITE(lunout,*)'Couplage Thermiques/Emanuel seulement si T<0'
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
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)
329	endif
330	RETURN
331	END
332

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

Download in other formats:

Original Format