Context Navigation

wake.F90 @ 4118

Last change on this file since 4118 was 4089, checked in by fhourdin, 3 years ago
Reecriture des thermiques
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: 68.1 KB

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

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

Download in other formats:

Original Format