Context Navigation

wake.F90 @ 4085

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

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

Download in other formats:

Original Format