Context Navigation

wake.F90 @ 4444

Last change on this file since 4444 was 4143, checked in by dcugnet, 3 years ago

Some variables are renamed or replaced by direct equivalents:
- iso_indnum -> tracers(:)%iso_iName
- niso_possibles -> niso
- iqiso -> iqIsoPha ; index_trac -> itZonIso
- ok_iso_verif -> isoCheck
- ntraceurs_zone -> nzone ; ntraciso -> ntiso
- qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio
Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times.
Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18).
The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).

File size: 81.6 KB

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

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

Download in other formats:

Original Format