Context Navigation

lmdz_wake.F90 @ 5454

Last change on this file since 5454 was 5285, checked in by abarral, 2 months ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 81.9 KB

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

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

Download in other formats:

Original Format