Context Navigation

wake.F90 @ 2473

Last change on this file since 2473 was 2467, checked in by idelkadi, 8 years ago
Utilisiation de la routine getin pour lire des parametres de la parametrization de la convection et nuages. Passage des fichiers conv_param.data et ep_param.data a conv_param.def et de wake_param.data a wake_param.def
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 53.1 KB

Line
1
2	! $Id: wake.F90 2467 2016-03-15 10:00:01Z fhourdin $
3
4	SUBROUTINE wake(p, ph, pi, dtime, sigd_con, &
5	te0, qe0, omgb, &
6	dtdwn, dqdwn, amdwn, amup, dta, dqa, &
7	wdtpbl, wdqpbl, udtpbl, udqpbl, &
8	deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, &
9	dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, &
10	dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
11	d_deltat_gw, d_deltatw2, d_deltaqw2)
12
13
14	! **************************************************************
15	! *
16	! WAKE *
17	! retour a un Pupper fixe *
18	! *
19	! written by : GRANDPEIX Jean-Yves 09/03/2000 *
20	! modified by : ROEHRIG Romain 01/29/2007 *
21	! **************************************************************
22
23	USE IOIPSL, ONLY : getin
24	USE dimphy
25	use mod_phys_lmdz_para
26	USE print_control_mod, ONLY: prt_level
27	IMPLICIT NONE
28	! ============================================================================
29
30
31	! But : Decrire le comportement des poches froides apparaissant dans les
32	! grands systemes convectifs, et fournir l'energie disponible pour
33	! le declenchement de nouvelles colonnes convectives.
34
35	! Variables d'etat : deltatw : ecart de temperature wake-undisturbed
36	! area
37	! deltaqw : ecart d'humidite wake-undisturbed area
38	! sigmaw : fraction d'aire occupee par la poche.
39
40	! Variable de sortie :
41
42	! wape : WAke Potential Energy
43	! fip : Front Incident Power (W/m2) - ALP
44	! gfl : Gust Front Length per unit area (m-1)
45	! dtls : large scale temperature tendency due to wake
46	! dqls : large scale humidity tendency due to wake
47	! hw : hauteur de la poche
48	! dp_omgb : vertical gradient of large scale omega
49	! wdens : densite de poches
50	! omgbdth: flux of Delta_Theta transported by LS omega
51	! dtKE : differential heating (wake - unpertubed)
52	! dqKE : differential moistening (wake - unpertubed)
53	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
54	! dp_deltomg : vertical gradient of omg (s-1)
55	! spread : spreading term in dt_wake and dq_wake
56	! deltatw : updated temperature difference (T_w-T_u).
57	! deltaqw : updated humidity difference (q_w-q_u).
58	! sigmaw : updated wake fractional area.
59	! d_deltat_gw : delta T tendency due to GW
60
61	! Variables d'entree :
62
63	! aire : aire de la maille
64	! te0 : temperature dans l'environnement (K)
65	! qe0 : humidite dans l'environnement (kg/kg)
66	! omgb : vitesse verticale moyenne sur la maille (Pa/s)
67	! dtdwn: source de chaleur due aux descentes (K/s)
68	! dqdwn: source d'humidite due aux descentes (kg/kg/s)
69	! dta : source de chaleur due courants satures et detrain (K/s)
70	! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
71	! amdwn: flux de masse total des descentes, par unite de
72	! surface de la maille (kg/m2/s)
73	! amup : flux de masse total des ascendances, par unite de
74	! surface de la maille (kg/m2/s)
75	! p : pressions aux milieux des couches (Pa)
76	! ph : pressions aux interfaces (Pa)
77	! pi : (p/p_0)**kapa (adim)
78	! dtime: increment temporel (s)
79
80	! Variables internes :
81
82	! rhow : masse volumique de la poche froide
83	! rho : environment density at P levels
84	! rhoh : environment density at Ph levels
85	! te : environment temperature \| may change within
86	! qe : environment humidity \| sub-time-stepping
87	! the : environment potential temperature
88	! thu : potential temperature in undisturbed area
89	! tu : temperature in undisturbed area
90	! qu : humidity in undisturbed area
91	! dp_omgb: vertical gradient og LS omega
92	! omgbw : wake average vertical omega
93	! dp_omgbw: vertical gradient of omgbw
94	! omgbdq : flux of Delta_q transported by LS omega
95	! dth : potential temperature diff. wake-undist.
96	! th1 : first pot. temp. for vertical advection (=thu)
97	! th2 : second pot. temp. for vertical advection (=thw)
98	! q1 : first humidity for vertical advection
99	! q2 : second humidity for vertical advection
100	! d_deltatw : terme de redistribution pour deltatw
101	! d_deltaqw : terme de redistribution pour deltaqw
102	! deltatw0 : deltatw initial
103	! deltaqw0 : deltaqw initial
104	! hw0 : hw initial
105	! sigmaw0: sigmaw initial
106	! amflux : horizontal mass flux through wake boundary
107	! wdens_ref: initial number of wakes per unit area (3D) or per
108	! unit length (2D), at the beginning of each time step
109	! Tgw : 1 sur la période de onde de gravité
110	! Cgw : vitesse de propagation de onde de gravité
111	! LL : distance entre 2 poches
112
113	! -------------------------------------------------------------------------
114	! Déclaration de variables
115	! -------------------------------------------------------------------------
116
117	include "YOMCST.h"
118	include "cvthermo.h"
119
120	! Arguments en entree
121	! --------------------
122
123	REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi
124	REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph, omgb
125	REAL, INTENT(IN) :: dtime
126	REAL, DIMENSION (klon, klev), INTENT(IN) :: te0, qe0
127	REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn
128	REAL, DIMENSION (klon, klev), INTENT(IN) :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED
129	REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup
130	REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa
131	REAL, DIMENSION (klon), INTENT(IN) :: sigd_con
132
133	!
134	! Input/Output
135	REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw
136	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw
137
138	! Sorties
139	! --------
140
141	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth
142	REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu
143	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls
144	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke
145	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtpbl, dqpbl
146	REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread
147	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2
148	REAL, DIMENSION (klon, klev+1), INTENT(OUT) :: omgbdth, omg
149	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg
150	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw
151	REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar
152	REAL, DIMENSION (klon), INTENT(OUT) :: wdens
153	INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw
154
155	! Variables internes
156	! -------------------
157
158	! Variables à fixer
159	REAL alon
160	LOGICAL, SAVE :: first = .TRUE.
161	!$OMP THREADPRIVATE(first)
162	REAL, SAVE :: stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol
163	!$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
164	REAL, SAVE :: stark_omp=0.33
165	REAL, SAVE :: alpk_omp=0.25
166	REAL, SAVE :: wdens_ref_omp=8.E-12
167	REAL, SAVE :: coefgw_omp=4.
168	!$OMP THREADPRIVATE(stark_omp,alpk_omp,wdens_ref_omp,coefgw_omp)
169
170	REAL delta_t_min
171	INTEGER nsub
172	REAL dtimesub
173	REAL sigmad, hwmin, wapecut
174	REAL :: sigmaw_max
175	REAL :: dens_rate
176	REAL wdens0
177	! IM 080208
178	LOGICAL, DIMENSION (klon) :: gwake
179
180	! Variables de sauvegarde
181	REAL, DIMENSION (klon, klev) :: deltatw0
182	REAL, DIMENSION (klon, klev) :: deltaqw0
183	REAL, DIMENSION (klon, klev) :: te, qe
184	REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
185
186	! Variables pour les GW
187	REAL, DIMENSION (klon) :: ll
188	REAL, DIMENSION (klon, klev) :: n2
189	REAL, DIMENSION (klon, klev) :: cgw
190	REAL, DIMENSION (klon, klev) :: tgw
191
192	! Variables liées au calcul de hw
193	REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
194	REAL, DIMENSION (klon) :: sum_dth
195	REAL, DIMENSION (klon) :: dthmin
196	REAL, DIMENSION (klon) :: z, dz, hw0
197	INTEGER, DIMENSION (klon) :: ktop, kupper
198
199	! Sub-timestep tendencies and related variables
200	REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
201	REAL d_te(klon, klev), d_qe(klon, klev)
202	REAL d_sigmaw(klon), alpha(klon)
203	REAL q0_min(klon), q1_min(klon)
204	LOGICAL wk_adv(klon), ok_qx_qw(klon)
205	REAL epsilon
206	DATA epsilon/1.E-15/
207
208	! Autres variables internes
209	INTEGER isubstep, k, i
210
211	REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
212	REAL, DIMENSION (klon) :: sum_dq, sum_rho
213	REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
214	REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
215	REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
216	REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
217
218	REAL, DIMENSION (klon, klev) :: rho, rhow
219	REAL, DIMENSION (klon, klev+1) :: rhoh
220	REAL, DIMENSION (klon, klev) :: rhow_moyen
221	REAL, DIMENSION (klon, klev) :: zh
222	REAL, DIMENSION (klon, klev+1) :: zhh
223	REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
224
225	REAL, DIMENSION (klon, klev) :: the, thu
226
227	! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
228
229	REAL, DIMENSION (klon, klev+1) :: omgbw
230	REAL, DIMENSION (klon) :: pupper
231	REAL, DIMENSION (klon) :: omgtop
232	REAL, DIMENSION (klon, klev) :: dp_omgbw
233	REAL, DIMENSION (klon) :: ztop, dztop
234	REAL, DIMENSION (klon, klev) :: alpha_up
235
236	REAL, DIMENSION (klon) :: rre1, rre2
237	REAL :: rrd1, rrd2
238	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
239	REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
240	REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
241	REAL, DIMENSION (klon, klev) :: omgbdq
242
243	REAL, DIMENSION (klon) :: ff, gg
244	REAL, DIMENSION (klon) :: wape2, cstar2, heff
245
246	REAL, DIMENSION (klon, klev) :: crep
247
248	REAL, DIMENSION (klon, klev) :: ppi
249
250	! cc nrlmd
251	REAL, DIMENSION (klon) :: death_rate, nat_rate
252	REAL, DIMENSION (klon, klev) :: entr
253	REAL, DIMENSION (klon, klev) :: detr
254
255	! -------------------------------------------------------------------------
256	! Initialisations
257	! -------------------------------------------------------------------------
258
259	! print*, 'wake initialisations'
260
261	! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
262	! -------------------------------------------------------------------------
263
264	DATA wapecut, sigmad, hwmin/5., .02, 10./
265	! cc nrlmd
266	DATA sigmaw_max/0.4/
267	DATA dens_rate/0.1/
268	! cc
269	! Longueur de maille (en m)
270	! -------------------------------------------------------------------------
271
272	! ALON = 3.e5
273	alon = 1.E6
274
275
276	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
277
278	! coefgw : Coefficient pour les ondes de gravité
279	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
280	! wdens : Densité de poche froide par maille
281	! -------------------------------------------------------------------------
282
283	! cc nrlmd coefgw=10
284	! coefgw=1
285	! wdens0 = 1.0/(alon**2)
286	! cc nrlmd wdens = 1.0/(alon**2)
287	! cc nrlmd stark = 0.50
288	! CRtest
289	! cc nrlmd alpk=0.1
290	! alpk = 1.0
291	! alpk = 0.5
292	! alpk = 0.05
293
294	if (first) then
295	crep_upper = 0.9
296	crep_sol = 1.0
297
298	! cc nrlmd Lecture du fichier wake_param.data
299	!$OMP MASTER
300	CALL getin('stark',stark_omp)
301	CALL getin('alpk',alpk_omp)
302	CALL getin('wdens_ref',wdens_ref_omp)
303	CALL getin('coefgw',coefgw_omp)
304	!$OMP END MASTER
305	!$OMP BARRIER
306	stark=stark_omp
307	alpk=alpk_omp
308	wdens_ref=wdens_ref_omp
309	coefgw=coefgw_omp
310	first=.false.
311	endif
312
313	! Initialisation de toutes des densites a wdens_ref.
314	! Les densites peuvent evoluer si les poches debordent
315	! (voir au tout debut de la boucle sur les substeps)
316	wdens = wdens_ref
317
318	! print*,'stark',stark
319	! print*,'alpk',alpk
320	! print*,'wdens',wdens
321	! print*,'coefgw',coefgw
322	! cc
323	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
324	! -------------------------------------------------------------------------
325
326	delta_t_min = 0.2
327
328	! 1. - Save initial values and initialize tendencies
329	! --------------------------------------------------
330
331	DO k = 1, klev
332	DO i = 1, klon
333	ppi(i, k) = pi(i, k)
334	deltatw0(i, k) = deltatw(i, k)
335	deltaqw0(i, k) = deltaqw(i, k)
336	te(i, k) = te0(i, k)
337	qe(i, k) = qe0(i, k)
338	dtls(i, k) = 0.
339	dqls(i, k) = 0.
340	d_deltat_gw(i, k) = 0.
341	d_te(i, k) = 0.
342	d_qe(i, k) = 0.
343	d_deltatw(i, k) = 0.
344	d_deltaqw(i, k) = 0.
345	! IM 060508 beg
346	d_deltatw2(i, k) = 0.
347	d_deltaqw2(i, k) = 0.
348	! IM 060508 end
349	END DO
350	END DO
351	! sigmaw1=sigmaw
352	! IF (sigd_con.GT.sigmaw1) THEN
353	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
354	! ENDIF
355	DO i = 1, klon
356	! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
357	sigmaw(i) = amax1(sigmaw(i), sigmad)
358	sigmaw(i) = amin1(sigmaw(i), 0.99)
359	sigmaw0(i) = sigmaw(i)
360	wape(i) = 0.
361	wape2(i) = 0.
362	d_sigmaw(i) = 0.
363	ktopw(i) = 0
364	END DO
365
366
367	! 2. - Prognostic part
368	! --------------------
369
370
371	! 2.1 - Undisturbed area and Wake integrals
372	! ---------------------------------------------------------
373
374	DO i = 1, klon
375	z(i) = 0.
376	ktop(i) = 0
377	kupper(i) = 0
378	sum_thu(i) = 0.
379	sum_tu(i) = 0.
380	sum_qu(i) = 0.
381	sum_thvu(i) = 0.
382	sum_dth(i) = 0.
383	sum_dq(i) = 0.
384	sum_rho(i) = 0.
385	sum_dtdwn(i) = 0.
386	sum_dqdwn(i) = 0.
387
388	av_thu(i) = 0.
389	av_tu(i) = 0.
390	av_qu(i) = 0.
391	av_thvu(i) = 0.
392	av_dth(i) = 0.
393	av_dq(i) = 0.
394	av_rho(i) = 0.
395	av_dtdwn(i) = 0.
396	av_dqdwn(i) = 0.
397	END DO
398
399	! Distance between wakes
400	DO i = 1, klon
401	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
402	END DO
403	! Potential temperatures and humidity
404	! ----------------------------------------------------------
405	DO k = 1, klev
406	DO i = 1, klon
407	! write(,)'wake 1',i,k,rd,te(i,k)
408	rho(i, k) = p(i, k)/(rd*te(i,k))
409	! write(,)'wake 2',rho(i,k)
410	IF (k==1) THEN
411	! write(,)'wake 3',i,k,rd,te(i,k)
412	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
413	! write(,)'wake 4',i,k,rd,te(i,k)
414	zhh(i, k) = 0
415	ELSE
416	! write(,)'wake 5',rd,(te(i,k)+te(i,k-1))
417	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
418	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
419	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
420	END IF
421	! write(,)'wake 7',ppi(i,k)
422	the(i, k) = te(i, k)/ppi(i, k)
423	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
424	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
425	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
426	! write(,)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
427	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
428	dth(i, k) = deltatw(i, k)/ppi(i, k)
429	END DO
430	END DO
431
432	DO k = 1, klev - 1
433	DO i = 1, klon
434	IF (k==1) THEN
435	n2(i, k) = 0
436	ELSE
437	n2(i, k) = amax1(0., -rg*2/the(i,k)rho(i,k)*(the(i,k+1)-the(i, &
438	k-1))/(p(i,k+1)-p(i,k-1)))
439	END IF
440	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
441
442	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
443	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
444	END DO
445	END DO
446
447	DO i = 1, klon
448	n2(i, klev) = 0
449	zh(i, klev) = 0
450	cgw(i, klev) = 0
451	tgw(i, klev) = 0
452	END DO
453
454	! Calcul de la masse volumique moyenne de la colonne (bdlmd)
455	! -----------------------------------------------------------------
456
457	DO k = 1, klev
458	DO i = 1, klon
459	epaisseur1(i, k) = 0.
460	epaisseur2(i, k) = 0.
461	END DO
462	END DO
463
464	DO i = 1, klon
465	epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
466	epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
467	rhow_moyen(i, 1) = rhow(i, 1)
468	END DO
469
470	DO k = 2, klev
471	DO i = 1, klon
472	epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
473	epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
474	rhow_moyen(i, k) = (rhow_moyen(i,k-1)epaisseur2(i,k-1)+rhow(i,k) &
475	epaisseur1(i,k))/epaisseur2(i, k)
476	END DO
477	END DO
478
479
480	! Choose an integration bound well above wake top
481	! -----------------------------------------------------------------
482
483	! Pupper = 50000. ! melting level
484	! Pupper = 60000.
485	! Pupper = 80000. ! essais pour case_e
486	DO i = 1, klon
487	pupper(i) = 0.6*ph(i, 1)
488	pupper(i) = max(pupper(i), 45000.)
489	! cc Pupper(i) = 60000.
490	END DO
491
492
493	! Determine Wake top pressure (Ptop) from buoyancy integral
494	! --------------------------------------------------------
495
496	! -1/ Pressure of the level where dth becomes less than delta_t_min.
497
498	DO i = 1, klon
499	ptop_provis(i) = ph(i, 1)
500	END DO
501	DO k = 2, klev
502	DO i = 1, klon
503
504	! IM v3JYG; ptop_provis(i).LT. ph(i,1)
505
506	IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
507	ptop_provis(i)==ph(i,1)) THEN
508	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
509	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
510	END IF
511	END DO
512	END DO
513
514	! -2/ dth integral
515
516	DO i = 1, klon
517	sum_dth(i) = 0.
518	dthmin(i) = -delta_t_min
519	z(i) = 0.
520	END DO
521
522	DO k = 1, klev
523	DO i = 1, klon
524	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
525	IF (dz(i)>0) THEN
526	z(i) = z(i) + dz(i)
527	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
528	dthmin(i) = amin1(dthmin(i), dth(i,k))
529	END IF
530	END DO
531	END DO
532
533	! -3/ height of triangle with area= sum_dth and base = dthmin
534
535	DO i = 1, klon
536	hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
537	hw0(i) = amax1(hwmin, hw0(i))
538	END DO
539
540	! -4/ now, get Ptop
541
542	DO i = 1, klon
543	z(i) = 0.
544	ptop(i) = ph(i, 1)
545	END DO
546
547	DO k = 1, klev
548	DO i = 1, klon
549	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
550	IF (dz(i)>0) THEN
551	z(i) = z(i) + dz(i)
552	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
553	END IF
554	END DO
555	END DO
556
557
558	! -5/ Determination de ktop et kupper
559
560	DO k = klev, 1, -1
561	DO i = 1, klon
562	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
563	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
564	END DO
565	END DO
566
567	! On evite kupper = 1 et kupper = klev
568	DO i = 1, klon
569	kupper(i) = max(kupper(i), 2)
570	kupper(i) = min(kupper(i), klev-1)
571	END DO
572
573
574	! -6/ Correct ktop and ptop
575
576	DO i = 1, klon
577	ptop_new(i) = ptop(i)
578	END DO
579	DO k = klev, 2, -1
580	DO i = 1, klon
581	IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
582	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
583	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
584	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
585	END IF
586	END DO
587	END DO
588
589	DO i = 1, klon
590	ptop(i) = ptop_new(i)
591	END DO
592
593	DO k = klev, 1, -1
594	DO i = 1, klon
595	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
596	END DO
597	END DO
598
599	! -5/ Set deltatw & deltaqw to 0 above kupper
600
601	DO k = 1, klev
602	DO i = 1, klon
603	IF (k>=kupper(i)) THEN
604	deltatw(i, k) = 0.
605	deltaqw(i, k) = 0.
606	END IF
607	END DO
608	END DO
609
610
611	! Vertical gradient of LS omega
612
613	DO k = 1, klev
614	DO i = 1, klon
615	IF (k<=kupper(i)) THEN
616	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
617	END IF
618	END DO
619	END DO
620
621	! Integrals (and wake top level number)
622	! --------------------------------------
623
624	! Initialize sum_thvu to 1st level virt. pot. temp.
625
626	DO i = 1, klon
627	z(i) = 1.
628	dz(i) = 1.
629	sum_thvu(i) = thu(i, 1)(1.+epsqu(i,1))*dz(i)
630	sum_dth(i) = 0.
631	END DO
632
633	DO k = 1, klev
634	DO i = 1, klon
635	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
636	IF (dz(i)>0) THEN
637	z(i) = z(i) + dz(i)
638	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
639	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
640	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
641	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsqu(i,k))*dz(i)
642	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
643	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
644	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
645	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
646	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
647	END IF
648	END DO
649	END DO
650
651	DO i = 1, klon
652	hw0(i) = z(i)
653	END DO
654
655
656	! 2.1 - WAPE and mean forcing computation
657	! ---------------------------------------
658
659	! ---------------------------------------
660
661	! Means
662
663	DO i = 1, klon
664	av_thu(i) = sum_thu(i)/hw0(i)
665	av_tu(i) = sum_tu(i)/hw0(i)
666	av_qu(i) = sum_qu(i)/hw0(i)
667	av_thvu(i) = sum_thvu(i)/hw0(i)
668	! av_thve = sum_thve/hw0
669	av_dth(i) = sum_dth(i)/hw0(i)
670	av_dq(i) = sum_dq(i)/hw0(i)
671	av_rho(i) = sum_rho(i)/hw0(i)
672	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
673	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
674
675	wape(i) = -rghw0(i)(av_dth(i)+eps(av_thu(i)av_dq(i)+av_dth(i)*av_qu(i &
676	)+av_dth(i)*av_dq(i)))/av_thvu(i)
677	END DO
678
679	! 2.2 Prognostic variable update
680	! ------------------------------
681
682	! Filter out bad wakes
683
684	DO k = 1, klev
685	DO i = 1, klon
686	IF (wape(i)<0.) THEN
687	deltatw(i, k) = 0.
688	deltaqw(i, k) = 0.
689	dth(i, k) = 0.
690	END IF
691	END DO
692	END DO
693
694	DO i = 1, klon
695	IF (wape(i)<0.) THEN
696	wape(i) = 0.
697	cstar(i) = 0.
698	hw(i) = hwmin
699	sigmaw(i) = amax1(sigmad, sigd_con(i))
700	fip(i) = 0.
701	gwake(i) = .FALSE.
702	ELSE
703	cstar(i) = starksqrt(2.wape(i))
704	gwake(i) = .TRUE.
705	END IF
706	END DO
707
708
709	! Check qx and qw positivity
710	! --------------------------
711	DO i = 1, klon
712	q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, &
713	1)+(1.-sigmaw(i))*deltaqw(i,1)))
714	END DO
715	DO k = 2, klev
716	DO i = 1, klon
717	q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, &
718	k)+(1.-sigmaw(i))*deltaqw(i,k)))
719	IF (q1_min(i)<=q0_min(i)) THEN
720	q0_min(i) = q1_min(i)
721	END IF
722	END DO
723	END DO
724
725	DO i = 1, klon
726	ok_qx_qw(i) = q0_min(i) >= 0.
727	alpha(i) = 1.
728	END DO
729
730	! C -----------------------------------------------------------------
731	! Sub-time-stepping
732	! -----------------
733
734	nsub = 10
735	dtimesub = dtime/nsub
736
737	! ------------------------------------------------------------
738	DO isubstep = 1, nsub
739	! ------------------------------------------------------------
740
741	! wk_adv is the logical flag enabling wake evolution in the time advance
742	! loop
743	DO i = 1, klon
744	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
745	END DO
746
747	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
748	! négatif de ktop à kupper --------
749	! cc On calcule pour cela une densité wdens0 pour laquelle on
750	! aurait un entrainement nul ---
751	DO i = 1, klon
752	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
753	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
754	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
755	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
756	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
757	wdens0 = (sigmaw(i)/(4.3.14))((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( &
758	ph(i,1)-pupper(i))cstar(i)))*(2)
759	IF (wdens(i)<=wdens0*1.1) THEN
760	wdens(i) = wdens0
761	END IF
762	! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
763	! c $ ,ph(i,1)-pupper(i)',
764	! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
765	! c $ ,ph(i,1)-pupper(i)
766	END IF
767	END DO
768
769	! cc nrlmd
770
771	DO i = 1, klon
772	IF (wk_adv(i)) THEN
773	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
774	sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
775	END IF
776	END DO
777	DO i = 1, klon
778	IF (wk_adv(i)) THEN
779	! cc nrlmd Introduction du taux de mortalité des poches et
780	! test sur sigmaw_max=0.4
781	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
782	IF (sigmaw(i)>=sigmaw_max) THEN
783	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
784	ELSE
785	death_rate(i) = 0.
786	END IF
787	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
788	dtimesub
789	! $ - nat_rate(i)sigmaw(i)dtimesub
790	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
791	! c $ death_rate(i),ktop(i),kupper(i)',
792	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
793	! c $ death_rate(i),ktop(i),kupper(i)
794
795	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
796	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
797	! wdens = wdens0/(10.*sigmaw)
798	! sigmaw =max(sigmaw,sigd_con)
799	! sigmaw =max(sigmaw,sigmad)
800	END IF
801	END DO
802
803
804	! calcul de la difference de vitesse verticale poche - zone non perturbee
805	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
806	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on
807	! definit
808	! IM 060208 au niveau k=1..?
809	DO k = 1, klev
810	DO i = 1, klon
811	IF (wk_adv(i)) THEN !!! nrlmd
812	dp_deltomg(i, k) = 0.
813	END IF
814	END DO
815	END DO
816	DO k = 1, klev + 1
817	DO i = 1, klon
818	IF (wk_adv(i)) THEN !!! nrlmd
819	omg(i, k) = 0.
820	END IF
821	END DO
822	END DO
823
824	DO i = 1, klon
825	IF (wk_adv(i)) THEN
826	z(i) = 0.
827	omg(i, 1) = 0.
828	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
829	END IF
830	END DO
831
832	DO k = 2, klev
833	DO i = 1, klon
834	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
835	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
836	z(i) = z(i) + dz(i)
837	dp_deltomg(i, k) = dp_deltomg(i, 1)
838	omg(i, k) = dp_deltomg(i, 1)*z(i)
839	END IF
840	END DO
841	END DO
842
843	DO i = 1, klon
844	IF (wk_adv(i)) THEN
845	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
846	ztop(i) = z(i) + dztop(i)
847	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
848	END IF
849	END DO
850
851	! -----------------
852	! From m/s to Pa/s
853	! -----------------
854
855	DO i = 1, klon
856	IF (wk_adv(i)) THEN
857	omgtop(i) = -rho(i, ktop(i))rgomgtop(i)
858	dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
859	END IF
860	END DO
861
862	DO k = 1, klev
863	DO i = 1, klon
864	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
865	omg(i, k) = -rho(i, k)rgomg(i, k)
866	dp_deltomg(i, k) = dp_deltomg(i, 1)
867	END IF
868	END DO
869	END DO
870
871	! raccordement lineaire de omg de ptop a pupper
872
873	DO i = 1, klon
874	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
875	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
876	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
877	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
878	(ptop(i)-pupper(i))
879	END IF
880	END DO
881
882	! c DO i=1,klon
883	! c print*,'Pente entre 0 et kupper (référence)'
884	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
885	! c print*,'Pente entre ktop et kupper'
886	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
887	! c ENDDO
888	! c
889	DO k = 1, klev
890	DO i = 1, klon
891	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
892	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
893	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
894	END IF
895	END DO
896	END DO
897	! cc nrlmd
898	! c DO i=1,klon
899	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
900	! c END DO
901	! cc
902
903
904	! -- Compute wake average vertical velocity omgbw
905
906
907	DO k = 1, klev + 1
908	DO i = 1, klon
909	IF (wk_adv(i)) THEN
910	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
911	END IF
912	END DO
913	END DO
914	! -- and its vertical gradient dp_omgbw
915
916	DO k = 1, klev
917	DO i = 1, klon
918	IF (wk_adv(i)) THEN
919	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
920	END IF
921	END DO
922	END DO
923
924	! -- Upstream coefficients for omgb velocity
925	! -- (alpha_up(k) is the coefficient of the value at level k)
926	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
927	DO k = 1, klev
928	DO i = 1, klon
929	IF (wk_adv(i)) THEN
930	alpha_up(i, k) = 0.
931	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
932	END IF
933	END DO
934	END DO
935
936	! Matrix expressing [The,deltatw] from [Th1,Th2]
937
938	DO i = 1, klon
939	IF (wk_adv(i)) THEN
940	rre1(i) = 1. - sigmaw(i)
941	rre2(i) = sigmaw(i)
942	END IF
943	END DO
944	rrd1 = -1.
945	rrd2 = 1.
946
947	! -- Get [Th1,Th2], dth and [q1,q2]
948
949	DO k = 1, klev
950	DO i = 1, klon
951	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
952	dth(i, k) = deltatw(i, k)/ppi(i, k)
953	th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
954	th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
955	q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
956	q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
957	END IF
958	END DO
959	END DO
960
961	DO i = 1, klon
962	IF (wk_adv(i)) THEN !!! nrlmd
963	d_th1(i, 1) = 0.
964	d_th2(i, 1) = 0.
965	d_dth(i, 1) = 0.
966	d_q1(i, 1) = 0.
967	d_q2(i, 1) = 0.
968	d_dq(i, 1) = 0.
969	END IF
970	END DO
971
972	DO k = 2, klev
973	DO i = 1, klon
974	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
975	d_th1(i, k) = th1(i, k-1) - th1(i, k)
976	d_th2(i, k) = th2(i, k-1) - th2(i, k)
977	d_dth(i, k) = dth(i, k-1) - dth(i, k)
978	d_q1(i, k) = q1(i, k-1) - q1(i, k)
979	d_q2(i, k) = q2(i, k-1) - q2(i, k)
980	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
981	END IF
982	END DO
983	END DO
984
985	DO i = 1, klon
986	IF (wk_adv(i)) THEN
987	omgbdth(i, 1) = 0.
988	omgbdq(i, 1) = 0.
989	END IF
990	END DO
991
992	DO k = 2, klev
993	DO i = 1, klon
994	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
995	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
996	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
997	END IF
998	END DO
999	END DO
1000
1001	! -----------------------------------------------------------------
1002	DO k = 1, klev
1003	DO i = 1, klon
1004	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1005	! -----------------------------------------------------------------
1006
1007	! Compute redistribution (advective) term
1008
1009	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1010	(rrd1omg(i,k)sigmaw(i)d_th1(i,k)-rrd2omg(i,k+1)*(1.-sigmaw( &
1011	i))d_th2(i,k+1)-(1.-alpha_up(i,k))omgbdth(i,k)-alpha_up(i,k+1)* &
1012	omgbdth(i,k+1))*ppi(i, k)
1013	! print*,'d_deltatw=',d_deltatw(i,k)
1014
1015	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1016	(rrd1omg(i,k)sigmaw(i)d_q1(i,k)-rrd2omg(i,k+1)*(1.-sigmaw( &
1017	i))d_q2(i,k+1)-(1.-alpha_up(i,k))omgbdq(i,k)-alpha_up(i,k+1)* &
1018	omgbdq(i,k+1))
1019	! print*,'d_deltaqw=',d_deltaqw(i,k)
1020
1021	! and increment large scale tendencies
1022
1023
1024
1025
1026	! C
1027	! -----------------------------------------------------------------
1028	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i, &
1029	k)-rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
1030	k+1)) & ! cc nrlmd $
1031	! -sigmaw(i)(1.-sigmaw(i))dth(i,k)*dp_deltomg(i,k)
1032	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
1033	k)-ph(i,k+1)) & ! cc
1034	)*ppi(i, k)
1035
1036	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i, &
1037	k)-rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
1038	k+1)) & ! cc nrlmd $
1039	! -sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*dp_deltomg(i,k)
1040	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i, &
1041	k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
1042	)
1043	! cc nrlmd
1044	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1045	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i, &
1046	k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k)
1047
1048	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i, &
1049	k)/(ph(i,k)-ph(i,k+1))))
1050
1051	END IF
1052	! cc
1053	END DO
1054	END DO
1055	! ------------------------------------------------------------------
1056
1057	! Increment state variables
1058
1059	DO k = 1, klev
1060	DO i = 1, klon
1061	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1062	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1063	! cc
1064
1065
1066
1067	! Coefficient de répartition
1068
1069	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1070	(ph(i,kupper(i))-ph(i,1))
1071	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i &
1072	,kupper(i)))
1073
1074
1075	! Reintroduce compensating subsidence term.
1076
1077	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1078	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1079	! . /(1-sigmaw)
1080	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1081	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1082	! . /(1-sigmaw)
1083
1084	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
1085	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
1086	! . /(1-sigmaw)
1087	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
1088	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
1089	! . /(1-sigmaw)
1090
1091	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1092	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1093	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1094
1095	!jyg<
1096	!!
1097	!!---------------------------------------------------------------
1098	!! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
1099	!! is accounted for by the PBL and the Thermals schemes. It is now set to zero
1100	!! within the Wake scheme.
1101	!!---------------------------------------------------------------
1102	dtPBL(i,k) = 0.
1103	dqPBL(i,k) = 0.
1104	!
1105	!! dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
1106	!! dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
1107	!
1108	!! dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
1109	!! dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
1110	! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
1111	!>jyg
1112	!
1113
1114	! cc nrlmd Prise en compte du taux de mortalité
1115	! cc Définitions de entr, detr
1116	detr(i, k) = 0.
1117
1118	entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1119	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1120
1121	spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1122	! cc spread(i,k) =
1123	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
1124	! cc $ sigmaw(i)
1125
1126
1127	! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
1128	! Jingmei
1129
1130	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1131	! & Tgw(i,k),deltatw(i,k)
1132	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
1133	dtimesub
1134	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1135	ff(i) = d_deltatw(i, k)/dtimesub
1136
1137	! Sans GW
1138
1139	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
1140
1141	! GW formule 1
1142
1143	! deltatw(k) = deltatw(k)+dtimesub*
1144	! $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
1145
1146	! GW formule 2
1147
1148	IF (dtimesub*tgw(i,k)<1.E-10) THEN
1149	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
1150	! $
1151	! -spread(i,k)*deltatw(i,k)
1152	-entr(i,k)deltatw(i,k)/sigmaw(i)-(death_rate(i)sigmaw( &
1153	i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1154	-tgw(i,k)*deltatw(i,k))
1155	ELSE
1156	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i, &
1157	k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $
1158	! -spread(i,k)*deltatw(i,k)
1159	-entr(i,k)deltatw(i,k)/sigmaw(i)-(death_rate(i)sigmaw( &
1160	i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
1161	-tgw(i,k)*deltatw(i,k))
1162	END IF
1163
1164	dth(i, k) = deltatw(i, k)/ppi(i, k)
1165
1166	gg(i) = d_deltaqw(i, k)/dtimesub
1167
1168	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $
1169	! -spread(i,k)*deltaqw(i,k))
1170	-entr(i,k)deltaqw(i,k)/sigmaw(i)-(death_rate(i)sigmaw(i)+detr( &
1171	i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
1172	! cc
1173
1174	! cc nrlmd
1175	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1176	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1177	! cc
1178	END IF
1179	END DO
1180	END DO
1181
1182
1183	! Scale tendencies so that water vapour remains positive in w and x.
1184
1185	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
1186	d_deltaqw, sigmaw, d_sigmaw, alpha)
1187
1188	! cc nrlmd
1189	! c print*,'alpha'
1190	! c do i=1,klon
1191	! c print*,alpha(i)
1192	! c end do
1193	! cc
1194	DO k = 1, klev
1195	DO i = 1, klon
1196	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1197	d_te(i, k) = alpha(i)*d_te(i, k)
1198	d_qe(i, k) = alpha(i)*d_qe(i, k)
1199	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1200	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1201	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1202	END IF
1203	END DO
1204	END DO
1205	DO i = 1, klon
1206	IF (wk_adv(i)) THEN
1207	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1208	END IF
1209	END DO
1210
1211	! Update large scale variables and wake variables
1212	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1213	! IM 060208 DO k = 1,kupper(i)
1214	DO k = 1, klev
1215	DO i = 1, klon
1216	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1217	dtls(i, k) = dtls(i, k) + d_te(i, k)
1218	dqls(i, k) = dqls(i, k) + d_qe(i, k)
1219	! cc nrlmd
1220	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1221	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1222	! cc
1223	END IF
1224	END DO
1225	END DO
1226	DO k = 1, klev
1227	DO i = 1, klon
1228	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1229	te(i, k) = te0(i, k) + dtls(i, k)
1230	qe(i, k) = qe0(i, k) + dqls(i, k)
1231	the(i, k) = te(i, k)/ppi(i, k)
1232	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1233	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1234	dth(i, k) = deltatw(i, k)/ppi(i, k)
1235	! c print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
1236	! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1237	END IF
1238	END DO
1239	END DO
1240	DO i = 1, klon
1241	IF (wk_adv(i)) THEN
1242	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1243	END IF
1244	END DO
1245
1246
1247	! Determine Ptop from buoyancy integral
1248	! ---------------------------------------
1249
1250	! - 1/ Pressure of the level where dth changes sign.
1251
1252	DO i = 1, klon
1253	IF (wk_adv(i)) THEN
1254	ptop_provis(i) = ph(i, 1)
1255	END IF
1256	END DO
1257
1258	DO k = 2, klev
1259	DO i = 1, klon
1260	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1261	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1262	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1263	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1264	END IF
1265	END DO
1266	END DO
1267
1268	! - 2/ dth integral
1269
1270	DO i = 1, klon
1271	IF (wk_adv(i)) THEN !!! nrlmd
1272	sum_dth(i) = 0.
1273	dthmin(i) = -delta_t_min
1274	z(i) = 0.
1275	END IF
1276	END DO
1277
1278	DO k = 1, klev
1279	DO i = 1, klon
1280	IF (wk_adv(i)) THEN
1281	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1282	IF (dz(i)>0) THEN
1283	z(i) = z(i) + dz(i)
1284	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1285	dthmin(i) = amin1(dthmin(i), dth(i,k))
1286	END IF
1287	END IF
1288	END DO
1289	END DO
1290
1291	! - 3/ height of triangle with area= sum_dth and base = dthmin
1292
1293	DO i = 1, klon
1294	IF (wk_adv(i)) THEN
1295	hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1296	hw(i) = amax1(hwmin, hw(i))
1297	END IF
1298	END DO
1299
1300	! - 4/ now, get Ptop
1301
1302	DO i = 1, klon
1303	IF (wk_adv(i)) THEN !!! nrlmd
1304	ktop(i) = 0
1305	z(i) = 0.
1306	END IF
1307	END DO
1308
1309	DO k = 1, klev
1310	DO i = 1, klon
1311	IF (wk_adv(i)) THEN
1312	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
1313	IF (dz(i)>0) THEN
1314	z(i) = z(i) + dz(i)
1315	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
1316	ktop(i) = k
1317	END IF
1318	END IF
1319	END DO
1320	END DO
1321
1322	! 4.5/Correct ktop and ptop
1323
1324	DO i = 1, klon
1325	IF (wk_adv(i)) THEN
1326	ptop_new(i) = ptop(i)
1327	END IF
1328	END DO
1329
1330	DO k = klev, 2, -1
1331	DO i = 1, klon
1332	! IM v3JYG; IF (k .GE. ktop(i)
1333	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1334	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1335	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
1336	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1337	END IF
1338	END DO
1339	END DO
1340
1341
1342	DO i = 1, klon
1343	IF (wk_adv(i)) THEN
1344	ptop(i) = ptop_new(i)
1345	END IF
1346	END DO
1347
1348	DO k = klev, 1, -1
1349	DO i = 1, klon
1350	IF (wk_adv(i)) THEN !!! nrlmd
1351	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
1352	END IF
1353	END DO
1354	END DO
1355
1356	! 5/ Set deltatw & deltaqw to 0 above kupper
1357
1358	DO k = 1, klev
1359	DO i = 1, klon
1360	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1361	deltatw(i, k) = 0.
1362	deltaqw(i, k) = 0.
1363	END IF
1364	END DO
1365	END DO
1366
1367
1368	! -------------Cstar computation---------------------------------
1369	DO i = 1, klon
1370	IF (wk_adv(i)) THEN !!! nrlmd
1371	sum_thu(i) = 0.
1372	sum_tu(i) = 0.
1373	sum_qu(i) = 0.
1374	sum_thvu(i) = 0.
1375	sum_dth(i) = 0.
1376	sum_dq(i) = 0.
1377	sum_rho(i) = 0.
1378	sum_dtdwn(i) = 0.
1379	sum_dqdwn(i) = 0.
1380
1381	av_thu(i) = 0.
1382	av_tu(i) = 0.
1383	av_qu(i) = 0.
1384	av_thvu(i) = 0.
1385	av_dth(i) = 0.
1386	av_dq(i) = 0.
1387	av_rho(i) = 0.
1388	av_dtdwn(i) = 0.
1389	av_dqdwn(i) = 0.
1390	END IF
1391	END DO
1392
1393	! Integrals (and wake top level number)
1394	! --------------------------------------
1395
1396	! Initialize sum_thvu to 1st level virt. pot. temp.
1397
1398	DO i = 1, klon
1399	IF (wk_adv(i)) THEN !!! nrlmd
1400	z(i) = 1.
1401	dz(i) = 1.
1402	sum_thvu(i) = thu(i, 1)(1.+epsqu(i,1))*dz(i)
1403	sum_dth(i) = 0.
1404	END IF
1405	END DO
1406
1407	DO k = 1, klev
1408	DO i = 1, klon
1409	IF (wk_adv(i)) THEN !!! nrlmd
1410	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1411	IF (dz(i)>0) THEN
1412	z(i) = z(i) + dz(i)
1413	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1414	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1415	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1416	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsqu(i,k))*dz(i)
1417	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1418	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1419	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1420	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1421	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1422	END IF
1423	END IF
1424	END DO
1425	END DO
1426
1427	DO i = 1, klon
1428	IF (wk_adv(i)) THEN !!! nrlmd
1429	hw0(i) = z(i)
1430	END IF
1431	END DO
1432
1433
1434	! - WAPE and mean forcing computation
1435	! ---------------------------------------
1436
1437	! ---------------------------------------
1438
1439	! Means
1440
1441	DO i = 1, klon
1442	IF (wk_adv(i)) THEN !!! nrlmd
1443	av_thu(i) = sum_thu(i)/hw0(i)
1444	av_tu(i) = sum_tu(i)/hw0(i)
1445	av_qu(i) = sum_qu(i)/hw0(i)
1446	av_thvu(i) = sum_thvu(i)/hw0(i)
1447	av_dth(i) = sum_dth(i)/hw0(i)
1448	av_dq(i) = sum_dq(i)/hw0(i)
1449	av_rho(i) = sum_rho(i)/hw0(i)
1450	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1451	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1452
1453	wape(i) = -rghw0(i)(av_dth(i)+eps(av_thu(i)av_dq(i)+av_dth(i)* &
1454	av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1455	END IF
1456	END DO
1457
1458	! Filter out bad wakes
1459
1460	DO k = 1, klev
1461	DO i = 1, klon
1462	IF (wk_adv(i)) THEN !!! nrlmd
1463	IF (wape(i)<0.) THEN
1464	deltatw(i, k) = 0.
1465	deltaqw(i, k) = 0.
1466	dth(i, k) = 0.
1467	END IF
1468	END IF
1469	END DO
1470	END DO
1471
1472	DO i = 1, klon
1473	IF (wk_adv(i)) THEN !!! nrlmd
1474	IF (wape(i)<0.) THEN
1475	wape(i) = 0.
1476	cstar(i) = 0.
1477	hw(i) = hwmin
1478	sigmaw(i) = max(sigmad, sigd_con(i))
1479	fip(i) = 0.
1480	gwake(i) = .FALSE.
1481	ELSE
1482	cstar(i) = starksqrt(2.wape(i))
1483	gwake(i) = .TRUE.
1484	END IF
1485	END IF
1486	END DO
1487
1488	END DO ! end sub-timestep loop
1489
1490	! -----------------------------------------------------------------
1491	! Get back to tendencies per second
1492
1493	DO k = 1, klev
1494	DO i = 1, klon
1495
1496	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1497	IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
1498	! cc
1499	dtls(i, k) = dtls(i, k)/dtime
1500	dqls(i, k) = dqls(i, k)/dtime
1501	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
1502	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
1503	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
1504	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
1505	! c $ ,death_rate(i)*sigmaw(i)
1506	END IF
1507	END DO
1508	END DO
1509
1510
1511	! ----------------------------------------------------------
1512	! Determine wake final state; recompute wape, cstar, ktop;
1513	! filter out bad wakes.
1514	! ----------------------------------------------------------
1515
1516	! 2.1 - Undisturbed area and Wake integrals
1517	! ---------------------------------------------------------
1518
1519	DO i = 1, klon
1520	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
1521	IF (ok_qx_qw(i)) THEN
1522	! cc
1523	z(i) = 0.
1524	sum_thu(i) = 0.
1525	sum_tu(i) = 0.
1526	sum_qu(i) = 0.
1527	sum_thvu(i) = 0.
1528	sum_dth(i) = 0.
1529	sum_dq(i) = 0.
1530	sum_rho(i) = 0.
1531	sum_dtdwn(i) = 0.
1532	sum_dqdwn(i) = 0.
1533
1534	av_thu(i) = 0.
1535	av_tu(i) = 0.
1536	av_qu(i) = 0.
1537	av_thvu(i) = 0.
1538	av_dth(i) = 0.
1539	av_dq(i) = 0.
1540	av_rho(i) = 0.
1541	av_dtdwn(i) = 0.
1542	av_dqdwn(i) = 0.
1543	END IF
1544	END DO
1545	! Potential temperatures and humidity
1546	! ----------------------------------------------------------
1547
1548	DO k = 1, klev
1549	DO i = 1, klon
1550	! cc nrlmd IF ( wk_adv(i)) THEN
1551	IF (ok_qx_qw(i)) THEN
1552	! cc
1553	rho(i, k) = p(i, k)/(rd*te(i,k))
1554	IF (k==1) THEN
1555	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
1556	zhh(i, k) = 0
1557	ELSE
1558	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
1559	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
1560	END IF
1561	the(i, k) = te(i, k)/ppi(i, k)
1562	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1563	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
1564	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
1565	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
1566	dth(i, k) = deltatw(i, k)/ppi(i, k)
1567	END IF
1568	END DO
1569	END DO
1570
1571	! Integrals (and wake top level number)
1572	! -----------------------------------------------------------
1573
1574	! Initialize sum_thvu to 1st level virt. pot. temp.
1575
1576	DO i = 1, klon
1577	! cc nrlmd IF ( wk_adv(i)) THEN
1578	IF (ok_qx_qw(i)) THEN
1579	! cc
1580	z(i) = 1.
1581	dz(i) = 1.
1582	sum_thvu(i) = thu(i, 1)(1.+epsqu(i,1))*dz(i)
1583	sum_dth(i) = 0.
1584	END IF
1585	END DO
1586
1587	DO k = 1, klev
1588	DO i = 1, klon
1589	! cc nrlmd IF ( wk_adv(i)) THEN
1590	IF (ok_qx_qw(i)) THEN
1591	! cc
1592	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1593	IF (dz(i)>0) THEN
1594	z(i) = z(i) + dz(i)
1595	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1596	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1597	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1598	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsqu(i,k))*dz(i)
1599	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1600	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1601	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1602	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1603	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1604	END IF
1605	END IF
1606	END DO
1607	END DO
1608
1609	DO i = 1, klon
1610	! cc nrlmd IF ( wk_adv(i)) THEN
1611	IF (ok_qx_qw(i)) THEN
1612	! cc
1613	hw0(i) = z(i)
1614	END IF
1615	END DO
1616
1617	! - WAPE and mean forcing computation
1618	! -------------------------------------------------------------
1619
1620	! Means
1621
1622	DO i = 1, klon
1623	! cc nrlmd IF ( wk_adv(i)) THEN
1624	IF (ok_qx_qw(i)) THEN
1625	! cc
1626	av_thu(i) = sum_thu(i)/hw0(i)
1627	av_tu(i) = sum_tu(i)/hw0(i)
1628	av_qu(i) = sum_qu(i)/hw0(i)
1629	av_thvu(i) = sum_thvu(i)/hw0(i)
1630	av_dth(i) = sum_dth(i)/hw0(i)
1631	av_dq(i) = sum_dq(i)/hw0(i)
1632	av_rho(i) = sum_rho(i)/hw0(i)
1633	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1634	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1635
1636	wape2(i) = -rghw0(i)(av_dth(i)+eps(av_thu(i)av_dq(i)+av_dth(i)* &
1637	av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
1638	END IF
1639	END DO
1640
1641	! Prognostic variable update
1642	! ------------------------------------------------------------
1643
1644	! Filter out bad wakes
1645
1646	DO k = 1, klev
1647	DO i = 1, klon
1648	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
1649	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
1650	! cc
1651	deltatw(i, k) = 0.
1652	deltaqw(i, k) = 0.
1653	dth(i, k) = 0.
1654	END IF
1655	END DO
1656	END DO
1657
1658
1659	DO i = 1, klon
1660	! cc nrlmd IF ( wk_adv(i)) THEN
1661	IF (ok_qx_qw(i)) THEN
1662	! cc
1663	IF (wape2(i)<0.) THEN
1664	wape2(i) = 0.
1665	cstar2(i) = 0.
1666	hw(i) = hwmin
1667	sigmaw(i) = amax1(sigmad, sigd_con(i))
1668	fip(i) = 0.
1669	gwake(i) = .FALSE.
1670	ELSE
1671	IF (prt_level>=10) PRINT *, 'wape2>0'
1672	cstar2(i) = starksqrt(2.wape2(i))
1673	gwake(i) = .TRUE.
1674	END IF
1675	END IF
1676	END DO
1677
1678	DO i = 1, klon
1679	! cc nrlmd IF ( wk_adv(i)) THEN
1680	IF (ok_qx_qw(i)) THEN
1681	! cc
1682	ktopw(i) = ktop(i)
1683	END IF
1684	END DO
1685
1686	DO i = 1, klon
1687	! cc nrlmd IF ( wk_adv(i)) THEN
1688	IF (ok_qx_qw(i)) THEN
1689	! cc
1690	IF (ktopw(i)>0 .AND. gwake(i)) THEN
1691
1692	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
1693	! cc heff = 600.
1694	! Utilisation de la hauteur hw
1695	! c heff = 0.7*hw
1696	heff(i) = hw(i)
1697
1698	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
1699	sqrt(sigmaw(i)wdens(i)3.14)
1700	fip(i) = alpk*fip(i)
1701	! jyg2
1702	ELSE
1703	fip(i) = 0.
1704	END IF
1705	END IF
1706	END DO
1707
1708	! Limitation de sigmaw
1709
1710	! cc nrlmd
1711	! DO i=1,klon
1712	! IF (OK_qx_qw(i)) THEN
1713	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
1714	! ENDIF
1715	! ENDDO
1716	! cc
1717	DO k = 1, klev
1718	DO i = 1, klon
1719
1720	! cc nrlmd On maintient désormais constant sigmaw en régime
1721	! permanent
1722	! cc IF ((sigmaw(i).GT.sigmaw_max).or.
1723	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1724	.NOT. ok_qx_qw(i)) THEN
1725	! cc
1726	dtls(i, k) = 0.
1727	dqls(i, k) = 0.
1728	deltatw(i, k) = 0.
1729	deltaqw(i, k) = 0.
1730	END IF
1731	END DO
1732	END DO
1733
1734	! cc nrlmd On maintient désormais constant sigmaw en régime permanent
1735	DO i = 1, klon
1736	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1737	.NOT. ok_qx_qw(i)) THEN
1738	wape(i) = 0.
1739	cstar(i) = 0.
1740	!!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
1741	!! hw(i) = hwmin !jyg
1742	!! sigmaw(i) = sigmad !jyg
1743	hw(i) = 0. !jyg
1744	sigmaw(i) = 0. !jyg
1745	fip(i) = 0.
1746	ELSE
1747	wape(i) = wape2(i)
1748	cstar(i) = cstar2(i)
1749	END IF
1750	! c print*,'wape wape2 ktopw OK_qx_qw =',
1751	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
1752	END DO
1753
1754
1755	RETURN
1756	END SUBROUTINE wake
1757
1758	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
1759	d_deltaqw, sigmaw, d_sigmaw, alpha)
1760	! ------------------------------------------------------
1761	! Dtermination du coefficient alpha tel que les tendances
1762	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
1763	! a une humidite positive dans la zone (x) et dans la zone (w).
1764	! ------------------------------------------------------
1765	IMPLICIT NONE
1766
1767	! Input
1768	REAL qe(nlon, nl), d_qe(nlon, nl)
1769	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
1770	REAL sigmaw(nlon), d_sigmaw(nlon)
1771	LOGICAL wk_adv(nlon)
1772	INTEGER nl, nlon
1773	! Output
1774	REAL alpha(nlon)
1775	! Internal variables
1776	REAL zeta(nlon, nl)
1777	REAL alpha1(nlon)
1778	REAL x, a, b, c, discrim
1779	REAL epsilon
1780	! DATA epsilon/1.e-15/
1781	INTEGER i,k
1782
1783	DO k = 1, nl
1784	DO i = 1, nlon
1785	IF (wk_adv(i)) THEN
1786	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
1787	zeta(i, k) = 0.
1788	ELSE
1789	zeta(i, k) = 1.
1790	END IF
1791	END IF
1792	END DO
1793	DO i = 1, nlon
1794	IF (wk_adv(i)) THEN
1795	x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
1796	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i)(deltaqw(i,k)+ &
1797	d_deltaqw(i,k))
1798	a = -d_sigmaw(i)*d_deltaqw(i, k)
1799	b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
1800	deltaqw(i, k)*d_sigmaw(i)
1801	c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
1802	discrim = bb - 4.a*c
1803	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
1804	IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
1805	alpha1(i) = 1.
1806	ELSE
1807	IF (x>=0.) THEN
1808	alpha1(i) = 1.
1809	ELSE
1810	IF (a>0.) THEN
1811	alpha1(i) = 0.9min((2.c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1812	))/(2.*a))
1813	ELSE IF (a==0.) THEN
1814	alpha1(i) = 0.9*(-c/b)
1815	ELSE
1816	! print*,'a,b,c discrim',a,b,c discrim
1817	alpha1(i) = 0.9max((2.c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
1818	))/(2.*a))
1819	END IF
1820	END IF
1821	END IF
1822	alpha(i) = min(alpha(i), alpha1(i))
1823	END IF
1824	END DO
1825	END DO
1826
1827	RETURN
1828	END SUBROUTINE wake_vec_modulation
1829
1830

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

Download in other formats:

Original Format