Context Navigation

wake.F90 @ 3208

Last change on this file since 3208 was 3208, checked in by jyg, 6 years ago
Implementation of a first crude model of the dynamic of wake population.
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: 71.4 KB

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

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

Download in other formats:

Original Format