Context Navigation

wake.F90 @ 2155

Last change on this file since 2155 was 2155, checked in by jyg, 10 years ago
wake.F90 : differential effect of PBL set to zero ; it is assumed to be accounted for by Turbulent diffusion and Thermal schemes
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: 52.0 KB

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

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

Download in other formats:

Original Format