Context Navigation

lmdz_wake.F90 @ 4825

Last change on this file since 4825 was 4825, checked in by fhourdin, 3 months ago
Correction plantage debug wake
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.4 KB

Line
1	MODULE lmdz_wake
2
3	! $Id: lmdz_wake.F90 4825 2024-02-16 10:56:24Z fhourdin $
4
5	CONTAINS
6
7	SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, &
8	tenv0, qe0, omgb, &
9	dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
10	sigd_con, Cin, &
11	deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, & ! state variables
12	dth, hw, wape, fip, gfl, &
13	dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
14	dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
15	d_deltat_gw, & ! tendencies
16	d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2) ! tendencies
17
18
19	! **************************************************************
20	! *
21	! WAKE *
22	! retour a un Pupper fixe *
23	! *
24	! written by : GRANDPEIX Jean-Yves 09/03/2000 *
25	! modified by : ROEHRIG Romain 01/29/2007 *
26	! **************************************************************
27
28
29	USE lmdz_wake_ini , ONLY : wake_ini
30	USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD
31	USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper
32	USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
33	USE lmdz_wake_ini , ONLY : ok_bug_gfl
34	USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold
35	USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
36	USE lmdz_wake_ini , ONLY : iflag_wk_profile
37	USE lmdz_wake_ini , ONLY : smallestreal
38
39
40	IMPLICIT NONE
41	! ============================================================================
42
43
44	! But : Decrire le comportement des poches froides apparaissant dans les
45	! grands systemes convectifs, et fournir l'energie disponible pour
46	! le declenchement de nouvelles colonnes convectives.
47
48	! State variables :
49	! deltatw : temperature difference between wake and off-wake regions
50	! deltaqw : specific humidity difference between wake and off-wake regions
51	! sigmaw : fractional area covered by wakes.
52	! asigmaw : fractional area covered by active wakes.
53	! wdens : number of wakes per unit area
54	! awdens : number of active wakes per unit area
55
56	! Variable de sortie :
57
58	! wape : WAke Potential Energy
59	! fip : Front Incident Power (W/m2) - ALP
60	! gfl : Gust Front Length per unit area (m-1)
61	! dtls : large scale temperature tendency due to wake
62	! dqls : large scale humidity tendency due to wake
63	! hw : wake top hight (given by hw*deltatw(1)/2=wape)
64	! dp_omgb : vertical gradient of large scale omega
65	! awdens : densite de poches actives
66	! wdens : densite de poches
67	! omgbdth: flux of Delta_Theta transported by LS omega
68	! dtKE : differential heating (wake - unpertubed)
69	! dqKE : differential moistening (wake - unpertubed)
70	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
71	! dp_deltomg : vertical gradient of omg (s-1)
72	! wkspread : spreading term in d_t_wake and d_q_wake
73	! deltatw : updated temperature difference (T_w-T_u).
74	! deltaqw : updated humidity difference (q_w-q_u).
75	! sigmaw : updated wake fractional area.
76	! asigmaw : updated active 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	! tenv0 : 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	! tenv : 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 periode de onde de gravite
130	! Cgw : vitesse de propagation de onde de gravite
131	! LL : distance entre 2 poches
132
133	! -------------------------------------------------------------------------
134	! Declaration de variables
135	! -------------------------------------------------------------------------
136
137
138	! Arguments en entree
139	! --------------------
140
141	INTEGER, INTENT(IN) :: klon,klev
142	INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf
143	REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi
144	REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph
145	REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb
146	REAL, INTENT(IN) :: dtime
147	REAL, DIMENSION (klon, klev), INTENT(IN) :: tenv0, qe0
148	REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn
149	REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup
150	REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa
151	REAL, DIMENSION (klon), INTENT(IN) :: wgen
152	REAL, DIMENSION (klon), INTENT(IN) :: sigd_con
153	REAL, DIMENSION (klon), INTENT(IN) :: Cin
154
155	!
156	! Input/Output
157	! State variables
158	REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw
159	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw
160	REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw
161	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens
162	REAL, DIMENSION (klon), INTENT(INOUT) :: awdens
163
164	! Sorties
165	! --------
166
167	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth
168	REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu
169	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls
170	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke
171	REAL, DIMENSION (klon, klev), INTENT(OUT) :: wkspread ! unused (jyg)
172	REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg
173	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg
174	REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar
175	INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw
176	! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
177	! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
178	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw
179	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2
180	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2
181
182	! Variables internes
183	! -------------------
184
185	! Variables a fixer
186
187	REAL :: delta_t_min
188	INTEGER :: nsub
189	REAL :: dtimesub
190	REAL :: wdens0
191	! IM 080208
192	LOGICAL, DIMENSION (klon) :: gwake
193
194	! Variables de sauvegarde
195	REAL, DIMENSION (klon, klev) :: deltatw0
196	REAL, DIMENSION (klon, klev) :: deltaqw0
197	REAL, DIMENSION (klon, klev) :: tenv, qe
198
199	! Variables liees a la dynamique de population 1
200	REAL, DIMENSION(klon) :: act
201	REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv
202	REAL, DIMENSION(klon) :: f_shear
203	REAL, DIMENSION(klon) :: drdt
204
205	! Variables liees a la dynamique de population 2
206	REAL, DIMENSION(klon) :: cont_fact
207
208	! Variables liees a la dynamique de population 3
209	REAL, DIMENSION(klon) :: arad_wk, irad_wk
210
211	!! REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col
212	REAL, DIMENSION(klon) :: wape1_act, wape2_act
213	LOGICAL, DIMENSION (klon) :: kill_wake
214	REAL :: drdt_pos
215	REAL :: tau_wk_inv_min
216	! Some components of the tendencies of state variables
217	REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
218	REAL, DIMENSION (klon) :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2
219	REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
220	REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
221
222	! Variables pour les GW
223	REAL, DIMENSION (klon) :: ll
224	REAL, DIMENSION (klon, klev) :: n2
225	REAL, DIMENSION (klon, klev) :: cgw
226	REAL, DIMENSION (klon, klev) :: tgw
227
228	! Variables liees au calcul de hw
229	REAL, DIMENSION (klon) :: ptop
230	REAL, DIMENSION (klon) :: sum_dth
231	REAL, DIMENSION (klon) :: dthmin
232	REAL, DIMENSION (klon) :: z, dz, hw0
233	INTEGER, DIMENSION (klon) :: ktop, kupper
234
235	! Variables liees au test de la forme triangulaire du profil de Delta_theta
236	REAL, DIMENSION (klon) :: sum_half_dth
237	REAL, DIMENSION (klon) :: dz_half
238
239	! Sub-timestep tendencies and related variables
240	REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw
241	REAL, DIMENSION (klon, klev) :: d_tenv, d_qe
242	REAL, DIMENSION (klon) :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
243	REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
244	REAL, DIMENSION (klon) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
245	REAL, DIMENSION (klon) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
246	REAL, DIMENSION (klon) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
247	REAL, DIMENSION (klon) :: agfl !! gust front length of active wakes
248	!! per unit area
249	REAL, DIMENSION (klon) :: alpha, alpha_tot
250	REAL, DIMENSION (klon) :: q0_min, q1_min
251	LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw
252
253	! Autres variables internes
254	INTEGER ::isubstep, k, i, igout
255
256	REAL :: wdensmin
257
258	REAL :: sigmaw_targ
259	REAL :: wdens_targ
260	REAL :: d_sigmaw_targ
261	REAL :: d_wdens_targ
262
263	REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
264	REAL, DIMENSION (klon) :: sum_dq, sum_rho
265	REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
266	REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
267	REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
268	REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
269
270	REAL, DIMENSION (klon, klev) :: rho, rhow
271	REAL, DIMENSION (klon, klev+1) :: rhoh
272	REAL, DIMENSION (klon, klev) :: rhow_moyen
273	REAL, DIMENSION (klon, klev) :: zh
274	REAL, DIMENSION (klon, klev+1) :: zhh
275	REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
276
277	REAL, DIMENSION (klon, klev) :: the, thu
278
279	REAL, DIMENSION (klon, klev) :: omgbw
280	REAL, DIMENSION (klon) :: pupper
281	REAL, DIMENSION (klon) :: omgtop
282	REAL, DIMENSION (klon, klev) :: dp_omgbw
283	REAL, DIMENSION (klon) :: ztop, dztop
284	REAL, DIMENSION (klon, klev) :: alpha_up
285
286	REAL, DIMENSION (klon) :: rre1, rre2
287	REAL :: rrd1, rrd2
288	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
289	REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
290	REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
291	REAL, DIMENSION (klon, klev) :: omgbdq
292
293	REAL, DIMENSION (klon) :: ff, gg
294	REAL, DIMENSION (klon) :: wape2, cstar2, heff
295
296	REAL, DIMENSION (klon, klev) :: crep
297
298	REAL, DIMENSION (klon, klev) :: ppi
299
300	! cc nrlmd
301	REAL, DIMENSION (klon) :: death_rate
302	!! REAL, DIMENSION (klon) :: nat_rate
303	REAL, DIMENSION (klon, klev) :: entr
304	REAL, DIMENSION (klon, klev) :: detr
305
306	REAL, DIMENSION(klon) :: sigmaw_in, asigmaw_in ! pour les prints
307	REAL, DIMENSION(klon) :: wdens_in, awdens_in ! pour les prints
308
309	!!! LOGICAL :: phys_sub=.true.
310	LOGICAL :: phys_sub=.false.
311
312	LOGICAL :: first_call=.true.
313
314
315	! -------------------------------------------------------------------------
316	! Initialisations
317	! -------------------------------------------------------------------------
318	! ALON = 3.e5
319	! alon = 1.E6
320
321	! Provisionnal; to be suppressed when f_shear is parameterized
322	f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear
323
324	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
325
326	! coefgw : Coefficient pour les ondes de gravite
327	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
328	! wdens : Densite surfacique de poche froide
329	! -------------------------------------------------------------------------
330
331	! cc nrlmd coefgw=10
332	! coefgw=1
333	! wdens0 = 1.0/(alon**2)
334	! cc nrlmd wdens = 1.0/(alon**2)
335	! cc nrlmd stark = 0.50
336	! CRtest
337	! cc nrlmd alpk=0.1
338	! alpk = 1.0
339	! alpk = 0.5
340	! alpk = 0.05
341	!
342	igout = klon/2+1/klon
343	!
344	! sub-time-stepping parameters
345	nsub = 10
346	dtimesub = dtime/nsub
347	!
348	IF (first_call) THEN
349	!!#define IOPHYS_WK
350	#undef IOPHYS_WK
351	#ifdef IOPHYS_WK
352	IF (phys_sub) THEN
353	call iophys_ini(dtimesub)
354	ELSE
355	call iophys_ini(dtime)
356	ENDIF
357	#endif
358	first_call = .false.
359	ENDIF !(first_call)
360
361	IF (iflag_wk_pop_dyn == 0) THEN
362	! Initialisation de toutes des densites a wdens_ref.
363	! Les densites peuvent evoluer si les poches debordent
364	! (voir au tout debut de la boucle sur les substeps)
365	!jyg<
366	!! wdens(:) = wdens_ref
367	DO i = 1,klon
368	wdens(i) = wdens_ref(znatsurf(i)+1)
369	ENDDO
370	!>jyg
371	ENDIF ! (iflag_wk_pop_dyn == 0)
372	!
373	IF (iflag_wk_pop_dyn >=1) THEN
374	IF (iflag_wk_pop_dyn == 3) THEN
375	wdensmin = wdensthreshold
376	ELSE
377	wdensmin = wdensinit
378	ENDIF
379	ENDIF
380
381	! print*,'stark',stark
382	! print*,'alpk',alpk
383	! print*,'wdens',wdens
384	! print*,'coefgw',coefgw
385	! cc
386	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
387	! -------------------------------------------------------------------------
388
389	delta_t_min = 0.2
390
391	! 1. - Save initial values, initialize tendencies, initialize output fields
392	! ------------------------------------------------------------------------
393
394	!jyg<
395	!! DO k = 1, klev
396	!! DO i = 1, klon
397	!! ppi(i, k) = pi(i, k)
398	!! deltatw0(i, k) = deltatw(i, k)
399	!! deltaqw0(i, k) = deltaqw(i, k)
400	!! tenv(i, k) = tenv0(i, k)
401	!! qe(i, k) = qe0(i, k)
402	!! dtls(i, k) = 0.
403	!! dqls(i, k) = 0.
404	!! d_deltat_gw(i, k) = 0.
405	!! d_tenv(i, k) = 0.
406	!! d_qe(i, k) = 0.
407	!! d_deltatw(i, k) = 0.
408	!! d_deltaqw(i, k) = 0.
409	!! ! IM 060508 beg
410	!! d_deltatw2(i, k) = 0.
411	!! d_deltaqw2(i, k) = 0.
412	!! ! IM 060508 end
413	!! END DO
414	!! END DO
415	ppi(:,:) = pi(:,:)
416	deltatw0(:,:) = deltatw(:,:)
417	deltaqw0(:,:) = deltaqw(:,:)
418	tenv(:,:) = tenv0(:,:)
419	qe(:,:) = qe0(:,:)
420	dtls(:,:) = 0.
421	dqls(:,:) = 0.
422	d_deltat_gw(:,:) = 0.
423	d_tenv(:,:) = 0.
424	d_qe(:,:) = 0.
425	d_deltatw(:,:) = 0.
426	d_deltaqw(:,:) = 0.
427	d_deltatw2(:,:) = 0.
428	d_deltaqw2(:,:) = 0.
429
430	d_sig_gen2(:) = 0.
431	d_sig_death2(:) = 0.
432	d_sig_col2(:) = 0.
433	d_sig_spread2(:)= 0.
434	d_asig_death2(:) = 0.
435	d_asig_iicol2(:) = 0.
436	d_asig_aicol2(:) = 0.
437	d_asig_spread2(:)= 0.
438	d_asig_bnd2(:) = 0.
439	d_asigmaw2(:) = 0.
440	!
441	d_dens_gen2(:) = 0.
442	d_dens_death2(:) = 0.
443	d_dens_col2(:) = 0.
444	d_dens_bnd2(:) = 0.
445	d_wdens2(:) = 0.
446	d_adens_bnd2(:) = 0.
447	d_awdens2(:) = 0.
448	d_adens_death2(:) = 0.
449	d_adens_icol2(:) = 0.
450	d_adens_acol2(:) = 0.
451
452	IF (iflag_wk_act == 0) THEN
453	act(:) = 0.
454	ELSEIF (iflag_wk_act == 1) THEN
455	act(:) = 1.
456	ENDIF
457
458	!! DO i = 1, klon
459	!! sigmaw_in(i) = sigmaw(i)
460	!! END DO
461	sigmaw_in(:) = sigmaw(:)
462	asigmaw_in(:) = asigmaw(:)
463	!>jyg
464	!
465	IF (iflag_wk_pop_dyn >= 1) THEN
466	awdens_in(:) = awdens(:)
467	wdens_in(:) = wdens(:)
468	!! wdens(:) = wdens(:) + wgen(:)*dtime
469	!! d_wdens2(:) = wgen(:)*dtime
470	!! ELSE
471	ENDIF ! (iflag_wk_pop_dyn >= 1)
472
473
474	! sigmaw1=sigmaw
475	! IF (sigd_con.GT.sigmaw1) THEN
476	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
477	! ENDIF
478	IF (iflag_wk_pop_dyn >= 1) THEN
479	DO i = 1, klon
480	d_dens_gen2(i) = 0.
481	d_dens_death2(i) = 0.
482	d_dens_col2(i) = 0.
483	d_awdens2(i) = 0.
484	IF (wdens(i) < wdensthreshold) THEN
485	!! wdens_targ = max(wdens(i),wdensmin)
486	wdens_targ = max(wdens(i),wdensinit)
487	d_dens_bnd2(i) = wdens_targ - wdens(i)
488	d_wdens2(i) = wdens_targ - wdens(i)
489	wdens(i) = wdens_targ
490	ELSE
491	d_dens_bnd2(i) = 0.
492	d_wdens2(i) = 0.
493	ENDIF !! (wdens(i) < wdensthreshold)
494	END DO
495	IF (iflag_wk_pop_dyn >= 2) THEN
496	DO i = 1, klon
497	IF (awdens(i) < wdensthreshold) THEN
498	!! wdens_targ = min(max(awdens(i),wdensmin),wdens(i))
499	wdens_targ = min(max(awdens(i),wdensinit),wdens(i))
500	d_adens_bnd2(i) = wdens_targ - awdens(i)
501	d_awdens2(i) = wdens_targ - awdens(i)
502	awdens(i) = wdens_targ
503	ELSE
504	wdens_targ = min(awdens(i), wdens(i))
505	d_adens_bnd2(i) = wdens_targ - awdens(i)
506	d_awdens2(i) = wdens_targ - awdens(i)
507	awdens(i) = wdens_targ
508	ENDIF
509	END DO
510	ENDIF ! (iflag_wk_pop_dyn >= 2)
511	ELSE
512	DO i = 1, klon
513	d_awdens2(i) = 0.
514	d_wdens2(i) = 0.
515	END DO
516	ENDIF ! (iflag_wk_pop_dyn >= 1)
517	!
518	DO i = 1, klon
519	sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
520	d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
521	d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
522	sigmaw(i) = sigmaw_targ
523	END DO
524	!
525	IF (iflag_wk_pop_dyn == 3) THEN
526	DO i = 1, klon
527	IF ((wdens(i)-awdens(i)) <= smallestreal) THEN
528	sigmaw_targ = sigmaw(i)
529	ELSE
530	sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
531	ENDIF
532	d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)
533	d_asigmaw2(i) = sigmaw_targ - asigmaw(i)
534	asigmaw(i) = sigmaw_targ
535	END DO
536	ENDIF ! (iflag_wk_pop_dyn == 3)
537
538	wape(:) = 0.
539	wape2(:) = 0.
540	d_sigmaw(:) = 0.
541	d_asigmaw(:) = 0.
542	ktopw(:) = 0
543	!
544	!<jyg
545	dth(:,:) = 0.
546	tu(:,:) = 0.
547	qu(:,:) = 0.
548	dtke(:,:) = 0.
549	dqke(:,:) = 0.
550	wkspread(:,:) = 0.
551	omgbdth(:,:) = 0.
552	omg(:,:) = 0.
553	dp_omgb(:,:) = 0.
554	dp_deltomg(:,:) = 0.
555	hw(:) = 0.
556	wape(:) = 0.
557	fip(:) = 0.
558	gfl(:) = 0.
559	cstar(:) = 0.
560	ktopw(:) = 0
561	!
562	! Vertical advection local variables
563	omgbw(:,:) = 0.
564	omgtop(:) = 0
565	dp_omgbw(:,:) = 0.
566	omgbdq(:,:) = 0.
567
568	!>jyg
569	!
570	IF (prt_level>=10) THEN
571	PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
572	PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
573	PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
574	PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
575	PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
576	PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
577	PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
578	PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
579	PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
580	ENDIF
581
582	! 2. - Prognostic part
583	! --------------------
584
585
586	! 2.1 - Undisturbed area and Wake integrals
587	! ---------------------------------------------------------
588
589	DO i = 1, klon
590	z(i) = 0.
591	ktop(i) = 0
592	kupper(i) = 0
593	sum_thu(i) = 0.
594	sum_tu(i) = 0.
595	sum_qu(i) = 0.
596	sum_thvu(i) = 0.
597	sum_dth(i) = 0.
598	sum_dq(i) = 0.
599	sum_rho(i) = 0.
600	sum_dtdwn(i) = 0.
601	sum_dqdwn(i) = 0.
602
603	av_thu(i) = 0.
604	av_tu(i) = 0.
605	av_qu(i) = 0.
606	av_thvu(i) = 0.
607	av_dth(i) = 0.
608	av_dq(i) = 0.
609	av_rho(i) = 0.
610	av_dtdwn(i) = 0.
611	av_dqdwn(i) = 0.
612	END DO
613
614	! Distance between wakes
615	DO i = 1, klon
616	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
617	END DO
618	! Potential temperatures and humidity
619	! ----------------------------------------------------------
620	DO k = 1, klev
621	DO i = 1, klon
622	! write(,)'wake 1',i,k,RD,tenv(i,k)
623	rho(i, k) = p(i, k)/(RD*tenv(i,k))
624	! write(,)'wake 2',rho(i,k)
625	IF (k==1) THEN
626	! write(,)'wake 3',i,k,rd,tenv(i,k)
627	rhoh(i, k) = ph(i, k)/(RD*tenv(i,k))
628	! write(,)'wake 4',i,k,rd,tenv(i,k)
629	zhh(i, k) = 0
630	ELSE
631	! write(,)'wake 5',rd,(tenv(i,k)+tenv(i,k-1))
632	rhoh(i, k) = ph(i, k)2./(RD(tenv(i,k)+tenv(i,k-1)))
633	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
634	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
635	END IF
636	! write(,)'wake 7',ppi(i,k)
637	the(i, k) = tenv(i, k)/ppi(i, k)
638	thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
639	tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i)
640	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
641	! write(,)'wake 8',(RD*(tenv(i,k)+deltatw(i,k)))
642	rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k)))
643	dth(i, k) = deltatw(i, k)/ppi(i, k)
644	END DO
645	END DO
646
647	DO k = 1, klev - 1
648	DO i = 1, klon
649	IF (k==1) THEN
650	n2(i, k) = 0
651	ELSE
652	n2(i, k) = amax1(0., -RG*2/the(i,k)rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
653	(p(i,k+1)-p(i,k-1)))
654	END IF
655	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
656
657	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
658	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
659	END DO
660	END DO
661
662	DO i = 1, klon
663	n2(i, klev) = 0
664	zh(i, klev) = 0
665	cgw(i, klev) = 0
666	tgw(i, klev) = 0
667	END DO
668
669	! Calcul de la masse volumique moyenne de la colonne (bdlmd)
670	! -----------------------------------------------------------------
671
672	DO k = 1, klev
673	DO i = 1, klon
674	epaisseur1(i, k) = 0.
675	epaisseur2(i, k) = 0.
676	END DO
677	END DO
678
679	DO i = 1, klon
680	epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.
681	epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.
682	rhow_moyen(i, 1) = rhow(i, 1)
683	END DO
684
685	DO k = 2, klev
686	DO i = 1, klon
687	epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG) + 1.
688	epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
689	rhow_moyen(i, k) = (rhow_moyen(i,k-1)epaisseur2(i,k-1)+rhow(i,k) &
690	epaisseur1(i,k))/epaisseur2(i, k)
691	END DO
692	END DO
693
694
695	! Choose an integration bound well above wake top
696	! -----------------------------------------------------------------
697
698	! Determine Wake top pressure (Ptop) from buoyancy integral
699	! --------------------------------------------------------
700	Do i=1, klon
701	wk_adv(i) = .True.
702	Enddo
703	Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
704	dth, hw0, rho, &
705	ktop, wk_adv)
706
707	IF (prt_level>=10) THEN
708	PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
709	ENDIF
710
711	! -5/ Set deltatw & deltaqw to 0 above kupper
712
713	DO k = 1, klev
714	DO i = 1, klon
715	IF (k>=kupper(i)) THEN
716	deltatw(i, k) = 0.
717	deltaqw(i, k) = 0.
718	d_deltatw2(i,k) = -deltatw0(i,k)
719	d_deltaqw2(i,k) = -deltaqw0(i,k)
720	END IF
721	END DO
722	END DO
723
724
725	! Vertical gradient of LS omega
726
727	DO k = 1, klev
728	DO i = 1, klon
729	IF (k<=kupper(i)) THEN
730	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
731	END IF
732	END DO
733	END DO
734
735	! Integrals (and wake top level number)
736	! --------------------------------------
737
738	! Initialize sum_thvu to 1st level virt. pot. temp.
739
740	DO i = 1, klon
741	z(i) = 1.
742	dz(i) = 1.
743	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
744	sum_dth(i) = 0.
745	END DO
746
747	DO k = 1, klev
748	DO i = 1, klon
749	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
750	IF (dz(i)>0) THEN
751	! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
752	z(i) = z(i) + dz(i)
753	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
754	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
755	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
756	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
757	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
758	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
759	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
760	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
761	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
762	END IF
763	END DO
764	END DO
765
766	DO i = 1, klon
767	hw0(i) = z(i)
768	END DO
769
770
771	! 2.1 - WAPE and mean forcing computation
772	! ---------------------------------------
773
774	! ---------------------------------------
775
776	! Means
777
778	DO i = 1, klon
779	av_thu(i) = sum_thu(i)/hw0(i)
780	av_tu(i) = sum_tu(i)/hw0(i)
781	av_qu(i) = sum_qu(i)/hw0(i)
782	av_thvu(i) = sum_thvu(i)/hw0(i)
783	! av_thve = sum_thve/hw0
784	av_dth(i) = sum_dth(i)/hw0(i)
785	av_dq(i) = sum_dq(i)/hw0(i)
786	av_rho(i) = sum_rho(i)/hw0(i)
787	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
788	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
789
790	wape(i) = -RGhw0(i)(av_dth(i)+ &
791	epsim1(av_thu(i)av_dq(i)+av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
792
793	END DO
794	#ifdef IOPHYS_WK
795	IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape)
796	#endif
797
798	! 2.2 Prognostic variable update
799	! ------------------------------
800
801	! Filter out bad wakes
802
803	DO k = 1, klev
804	DO i = 1, klon
805	IF (wape(i)<0.) THEN
806	deltatw(i, k) = 0.
807	deltaqw(i, k) = 0.
808	dth(i, k) = 0.
809	d_deltatw2(i,k) = -deltatw0(i,k)
810	d_deltaqw2(i,k) = -deltaqw0(i,k)
811	END IF
812	END DO
813	END DO
814
815	DO i = 1, klon
816	IF (wape(i)<0.) THEN
817	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
818	sigmaw_targ = max(sigmad, sigd_con(i))
819	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
820	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
821	sigmaw(i) = sigmaw_targ
822	ENDIF !! (wape(i)<0.)
823	ENDDO
824	!
825	IF (iflag_wk_pop_dyn == 3) THEN
826	DO i = 1, klon
827	IF (wape(i)<0.) THEN
828	sigmaw_targ = max(sigmad, sigd_con(i))
829	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
830	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
831	asigmaw(i) = sigmaw_targ
832	ENDIF !! (wape(i)<0.)
833	ENDDO
834	ENDIF !! (iflag_wk_pop_dyn == 3)
835
836	DO i = 1, klon
837	IF (wape(i)<0.) THEN
838	wape(i) = 0.
839	cstar(i) = 0.
840	hw(i) = hwmin
841	fip(i) = 0.
842	gwake(i) = .FALSE.
843	ELSE
844	hw(i) = hw0(i)
845	cstar(i) = starksqrt(2.wape(i))
846	gwake(i) = .TRUE.
847	END IF
848	END DO
849	!
850
851	! Check qx and qw positivity
852	! --------------------------
853	DO i = 1, klon
854	q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), &
855	(qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
856	END DO
857	DO k = 2, klev
858	DO i = 1, klon
859	q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
860	(qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
861	IF (q1_min(i)<=q0_min(i)) THEN
862	q0_min(i) = q1_min(i)
863	END IF
864	END DO
865	END DO
866
867	DO i = 1, klon
868	ok_qx_qw(i) = q0_min(i) >= 0.
869	alpha(i) = 1.
870	alpha_tot(i) = 1.
871	END DO
872
873	IF (prt_level>=10) THEN
874	PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
875	sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
876	ENDIF
877
878
879	! C -----------------------------------------------------------------
880	! Sub-time-stepping
881	! -----------------
882
883	! nsub and dtimesub definitions moved to begining of routine.
884	!! nsub = 10
885	!! dtimesub = dtime/nsub
886
887
888	! ------------------------------------------------------------------------
889	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
890	! ------------------------------------------------------------------------
891	!
892	wk_adv = .True.
893	DO isubstep = 1, nsub
894	!
895	! ------------------------------------------------------------------------
896	! wk_adv is the logical flag enabling wake evolution in the time advance
897	! loop
898
899	DO i = 1, klon
900	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
901	END DO
902	IF (prt_level>=10) THEN
903	PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
904	isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
905
906	ENDIF
907
908	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
909	! negatif de ktop a kupper --------
910	! cc On calcule pour cela une densite wdens0 pour laquelle on
911	! aurait un entrainement nul ---
912	!jyg<
913	! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
914	! les poches sont insuffisantes pour accueillir tout le flux de masse
915	! des descentes unsaturees. Nous faisons alors l'hypothese que la
916	! convection profonde cree directement de nouvelles poches, sans passer
917	! par les thermiques. La nouvelle valeur de wdens est alors imposee.
918
919	DO i = 1, klon
920	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
921	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
922	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
923	IF ( iflag_wk_profile == 0 ) THEN
924	omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
925	RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
926	ELSE
927	omg(i, kupper(i)+1)=0.
928	ENDIF
929	wdens0 = (sigmaw(i)/(4.3.14)) &
930	((1.-sigmaw(i))omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))cstar(i)))**(2)
931	IF (prt_level >= 10) THEN
932	print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
933	omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
934	ENDIF
935	IF (wdens(i)<=wdens0*1.1) THEN
936	IF (iflag_wk_pop_dyn >= 1) THEN
937	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
938	d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
939	ENDIF
940	wdens(i) = wdens0
941	END IF
942	END IF
943	END DO
944
945	IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN
946	!!--------------------------------------------------------
947	!!Bug : computing gfl and rad_wk before changing sigmaw
948	!! This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk
949	!! are computed within wake_popdyn
950	!!--------------------------------------------------------
951	DO i = 1, klon
952	IF (wk_adv(i)) THEN
953	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
954	rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
955	END IF
956	END DO
957	ENDIF ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl)
958	!!--------------------------------------------------------
959
960	DO i = 1, klon
961	IF (wk_adv(i)) THEN
962	sigmaw_targ = min(sigmaw(i), sigmaw_max)
963	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
964	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
965	sigmaw(i) = sigmaw_targ
966	END IF
967	END DO
968
969	IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN
970	!!--------------------------------------------------------
971	!!Fix : computing gfl and rad_wk after changing sigmaw
972	!!--------------------------------------------------------
973	DO i = 1, klon
974	IF (wk_adv(i)) THEN
975	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
976	rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
977	END IF
978	END DO
979	ENDIF ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl)
980	!!--------------------------------------------------------
981
982	IF (iflag_wk_pop_dyn >= 1) THEN
983	! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
984	! Here, it has to be set to zero.
985	death_rate(:) = 0.
986	ENDIF
987
988	IF (iflag_wk_pop_dyn >= 3) THEN
989	DO i = 1, klon
990	IF (wk_adv(i)) THEN
991	sigmaw_targ = min(asigmaw(i), sigmaw_max)
992	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
993	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
994	asigmaw(i) = sigmaw_targ
995	ENDIF
996	ENDDO
997	ENDIF
998
999	!!--------------------------------------------------------
1000	!!--------------------------------------------------------
1001	IF (iflag_wk_pop_dyn == 1) THEN
1002	!
1003	CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
1004	wdensmin, &
1005	dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
1006	d_awdens, d_wdens, d_sigmaw, &
1007	iflag_wk_act, wk_adv, cin, wape, &
1008	drdt, &
1009	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
1010	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
1011	d_wdens_targ, d_sigmaw_targ)
1012
1013
1014	!!--------------------------------------------------------
1015	ELSEIF (iflag_wk_pop_dyn == 2) THEN
1016	!
1017	CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
1018	wdensmin, &
1019	sigmaw, wdens, awdens, & !! state variables
1020	gfl, cstar, cin, wape, rad_wk, &
1021	d_sigmaw, d_wdens, d_awdens, & !! tendencies
1022	cont_fact, &
1023	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
1024	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
1025	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
1026	sigmaw=sigmaw-d_sigmaw
1027	wdens=wdens-d_wdens
1028	awdens=awdens-d_awdens
1029
1030	!!--------------------------------------------------------
1031	ELSEIF (iflag_wk_pop_dyn == 3) THEN
1032	#ifdef IOPHYS_WK
1033	IF (phys_sub) THEN
1034	CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop)
1035	CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw)
1036	CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw)
1037	CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens)
1038	CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens)
1039	CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk)
1040	CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk)
1041	CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk)
1042	ENDIF
1043	#endif
1044	!
1045	CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
1046	wdensmin, &
1047	sigmaw, asigmaw, wdens, awdens, & !! state variables
1048	gfl, agfl, cstar, cin, wape, &
1049	rad_wk, arad_wk, irad_wk, &
1050	d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies
1051	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
1052	d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
1053	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
1054	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
1055	sigmaw=sigmaw-d_sigmaw
1056	asigmaw=asigmaw-d_asigmaw
1057	wdens=wdens-d_wdens
1058	awdens=awdens-d_awdens
1059
1060	!!--------------------------------------------------------
1061	ELSEIF (iflag_wk_pop_dyn == 0) THEN
1062
1063	! cc nrlmd
1064
1065	DO i = 1, klon
1066	IF (wk_adv(i)) THEN
1067
1068	! cc nrlmd Introduction du taux de mortalite des poches et
1069	! test sur sigmaw_max=0.4
1070	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
1071	IF (sigmaw(i)>=sigmaw_max) THEN
1072	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
1073	ELSE
1074	death_rate(i) = 0.
1075	END IF
1076
1077	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
1078	dtimesub
1079	! $ - nat_rate(i)sigmaw(i)dtimesub
1080	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1081	! c $ death_rate(i),ktop(i),kupper(i)',
1082	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1083	! c $ death_rate(i),ktop(i),kupper(i)
1084
1085	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
1086	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
1087	! wdens = wdens0/(10.*sigmaw)
1088	! sigmaw =max(sigmaw,sigd_con)
1089	! sigmaw =max(sigmaw,sigmad)
1090	END IF
1091	END DO
1092
1093	ENDIF ! (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0)
1094	!!--------------------------------------------------------
1095	!!--------------------------------------------------------
1096
1097	#ifdef IOPHYS_WK
1098	IF (phys_sub) THEN
1099	CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens)
1100	CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens)
1101	CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw)
1102	CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw)
1103	ENDIF
1104	#endif
1105	! calcul de la difference de vitesse verticale poche - zone non perturbee
1106	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
1107	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
1108	! IM 060208 au niveau k=1...
1109	!JYG 161013 Correction : maintenant omg est dimensionne a klev.
1110	DO k = 1, klev
1111	DO i = 1, klon
1112	IF (wk_adv(i)) THEN !!! nrlmd
1113	dp_deltomg(i, k) = 0.
1114	END IF
1115	END DO
1116	END DO
1117	DO k = 1, klev
1118	DO i = 1, klon
1119	IF (wk_adv(i)) THEN !!! nrlmd
1120	omg(i, k) = 0.
1121	END IF
1122	END DO
1123	END DO
1124
1125	DO i = 1, klon
1126	IF (wk_adv(i)) THEN
1127	z(i) = 0.
1128	omg(i, 1) = 0.
1129	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
1130	END IF
1131	END DO
1132
1133	DO k = 2, klev
1134	DO i = 1, klon
1135	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1136	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
1137	z(i) = z(i) + dz(i)
1138	dp_deltomg(i, k) = dp_deltomg(i, 1)
1139	omg(i, k) = dp_deltomg(i, 1)*z(i)
1140	END IF
1141	END DO
1142	END DO
1143
1144	DO i = 1, klon
1145	IF (wk_adv(i)) THEN
1146	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
1147	ztop(i) = z(i) + dztop(i)
1148	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
1149	END IF
1150	END DO
1151
1152	IF (prt_level>=10) THEN
1153	PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
1154	PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
1155	omgtop(igout), ptop(igout), ktop(igout)
1156	ENDIF
1157
1158	! -----------------
1159	! From m/s to Pa/s
1160	! -----------------
1161
1162	DO i = 1, klon
1163	IF (wk_adv(i)) THEN
1164	omgtop(i) = -rho(i, ktop(i))RGomgtop(i)
1165	dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
1166	END IF
1167	END DO
1168
1169	DO k = 1, klev
1170	DO i = 1, klon
1171	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1172	omg(i, k) = -rho(i, k)RGomg(i, k)
1173	dp_deltomg(i, k) = dp_deltomg(i, 1)
1174	END IF
1175	END DO
1176	END DO
1177
1178	! raccordement lineaire de omg de ptop a pupper
1179
1180	DO i = 1, klon
1181	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
1182	IF ( iflag_wk_profile == 0 ) THEN
1183	omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
1184	RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
1185	ELSE
1186	omg(i, kupper(i)+1) = 0.
1187	ENDIF
1188	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
1189	(ptop(i)-pupper(i))
1190	END IF
1191	END DO
1192
1193	! c DO i=1,klon
1194	! c print*,'Pente entre 0 et kupper (reference)'
1195	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
1196	! c print*,'Pente entre ktop et kupper'
1197	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
1198	! c ENDDO
1199	! c
1200	DO k = 1, klev
1201	DO i = 1, klon
1202	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
1203	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
1204	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
1205	END IF
1206	END DO
1207	END DO
1208	!! print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
1209	! cc nrlmd
1210	! c DO i=1,klon
1211	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
1212	! c END DO
1213	! cc
1214
1215
1216	! -- Compute wake average vertical velocity omgbw
1217
1218
1219	DO k = 1, klev
1220	DO i = 1, klon
1221	IF (wk_adv(i)) THEN
1222	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
1223	END IF
1224	END DO
1225	END DO
1226	! -- and its vertical gradient dp_omgbw
1227
1228	DO k = 1, klev-1
1229	DO i = 1, klon
1230	IF (wk_adv(i)) THEN
1231	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
1232	END IF
1233	END DO
1234	END DO
1235	DO i = 1, klon
1236	IF (wk_adv(i)) THEN
1237	dp_omgbw(i, klev) = 0.
1238	END IF
1239	END DO
1240
1241	! -- Upstream coefficients for omgb velocity
1242	! -- (alpha_up(k) is the coefficient of the value at level k)
1243	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
1244	DO k = 1, klev
1245	DO i = 1, klon
1246	IF (wk_adv(i)) THEN
1247	alpha_up(i, k) = 0.
1248	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
1249	END IF
1250	END DO
1251	END DO
1252
1253	! Matrix expressing [The,deltatw] from [Th1,Th2]
1254
1255	DO i = 1, klon
1256	IF (wk_adv(i)) THEN
1257	rre1(i) = 1. - sigmaw(i)
1258	rre2(i) = sigmaw(i)
1259	END IF
1260	END DO
1261	rrd1 = -1.
1262	rrd2 = 1.
1263
1264	! -- Get [Th1,Th2], dth and [q1,q2]
1265
1266	DO k = 1, klev
1267	DO i = 1, klon
1268	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1269	dth(i, k) = deltatw(i, k)/ppi(i, k)
1270	th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
1271	th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
1272	q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
1273	q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
1274	END IF
1275	END DO
1276	END DO
1277
1278	DO i = 1, klon
1279	IF (wk_adv(i)) THEN !!! nrlmd
1280	d_th1(i, 1) = 0.
1281	d_th2(i, 1) = 0.
1282	d_dth(i, 1) = 0.
1283	d_q1(i, 1) = 0.
1284	d_q2(i, 1) = 0.
1285	d_dq(i, 1) = 0.
1286	END IF
1287	END DO
1288
1289	DO k = 2, klev
1290	DO i = 1, klon
1291	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1292	d_th1(i, k) = th1(i, k-1) - th1(i, k)
1293	d_th2(i, k) = th2(i, k-1) - th2(i, k)
1294	d_dth(i, k) = dth(i, k-1) - dth(i, k)
1295	d_q1(i, k) = q1(i, k-1) - q1(i, k)
1296	d_q2(i, k) = q2(i, k-1) - q2(i, k)
1297	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
1298	END IF
1299	END DO
1300	END DO
1301
1302	DO i = 1, klon
1303	IF (wk_adv(i)) THEN
1304	omgbdth(i, 1) = 0.
1305	omgbdq(i, 1) = 0.
1306	END IF
1307	END DO
1308
1309	DO k = 2, klev
1310	DO i = 1, klon
1311	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
1312	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
1313	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
1314	END IF
1315	END DO
1316	END DO
1317
1318	!! IF (prt_level>=10) THEN
1319	IF (prt_level>=10 .and. wk_adv(igout)) THEN
1320	PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
1321	PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
1322	PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
1323	PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
1324	ENDIF
1325
1326	! -----------------------------------------------------------------
1327	DO k = 1, klev-1
1328	DO i = 1, klon
1329	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1330	! -----------------------------------------------------------------
1331
1332	! Compute redistribution (advective) term
1333
1334	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1335	(rrd1omg(i,k)sigmaw(i)*d_th1(i,k) - &
1336	rrd2omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1)- &
1337	(1.-alpha_up(i,k))*omgbdth(i,k)- &
1338	alpha_up(i,k+1)omgbdth(i,k+1))ppi(i, k)
1339	! print*,'d_deltatw=', k, d_deltatw(i,k)
1340
1341	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1342	(rrd1omg(i,k)sigmaw(i)*d_q1(i,k)- &
1343	rrd2omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1)- &
1344	(1.-alpha_up(i,k))*omgbdq(i,k)- &
1345	alpha_up(i,k+1)*omgbdq(i,k+1))
1346	! print*,'d_deltaqw=', k, d_deltaqw(i,k)
1347
1348	! and increment large scale tendencies
1349
1350
1351
1352
1353	! C
1354	! -----------------------------------------------------------------
1355	d_tenv(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)- &
1356	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/ &
1357	(ph(i,k)-ph(i,k+1)) &
1358	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
1359	(ph(i,k)-ph(i,k+1)) )*ppi(i, k)
1360
1361	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)- &
1362	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/ &
1363	(ph(i,k)-ph(i,k+1)) &
1364	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
1365	(ph(i,k)-ph(i,k+1)) )
1366	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1367	d_tenv(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
1368
1369	d_qe(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
1370
1371	END IF
1372	! cc
1373	END DO
1374	END DO
1375	! ------------------------------------------------------------------
1376
1377	IF (prt_level>=10) THEN
1378	PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
1379	PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
1380	ENDIF
1381
1382	! Increment state variables
1383	!jyg<
1384	IF (iflag_wk_pop_dyn >= 1) THEN
1385	DO k = 1, klev
1386	DO i = 1, klon
1387	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1388	detr(i,k) = - d_sig_death(i) - d_sig_col(i)
1389	entr(i,k) = d_sig_gen(i)
1390	ENDIF
1391	ENDDO
1392	ENDDO
1393	ELSE ! (iflag_wk_pop_dyn >= 1)
1394	DO k = 1, klev
1395	DO i = 1, klon
1396	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1397	detr(i, k) = 0.
1398
1399	entr(i, k) = 0.
1400	ENDIF
1401	ENDDO
1402	ENDDO
1403	ENDIF ! (iflag_wk_pop_dyn >= 1)
1404
1405
1406
1407	DO k = 1, klev
1408	DO i = 1, klon
1409	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1410	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1411	! cc
1412
1413
1414
1415	! Coefficient de repartition
1416
1417	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1418	(ph(i,kupper(i))-ph(i,1))
1419	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
1420	(ph(i,1)-ph(i,kupper(i)))
1421
1422
1423	! Reintroduce compensating subsidence term.
1424
1425	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1426	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1427	! . /(1-sigmaw)
1428	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1429	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1430	! . /(1-sigmaw)
1431
1432	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
1433	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
1434	! . /(1-sigmaw)
1435	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
1436	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
1437	! . /(1-sigmaw)
1438
1439	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1440	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1441	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1442
1443	!
1444
1445	! cc nrlmd Prise en compte du taux de mortalite
1446	! cc Definitions de entr, detr
1447	!jyg<
1448	!! detr(i, k) = 0.
1449	!!
1450	!! entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1451	!! sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1452	!!
1453	entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
1454	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1455	!>jyg
1456	wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1457
1458	! cc wkspread(i,k) =
1459	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
1460	! cc $ sigmaw(i)
1461
1462
1463	! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
1464	! Jingmei
1465
1466	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1467	! & Tgw(i,k),deltatw(i,k)
1468	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
1469	dtimesub
1470	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1471	ff(i) = d_deltatw(i, k)/dtimesub
1472
1473	! Sans GW
1474
1475	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-wkspread(k)deltatw(k))
1476
1477	! GW formule 1
1478
1479	! deltatw(k) = deltatw(k)+dtimesub*
1480	! $ (ff+dtKE(k) - wkspread(k)deltatw(k)-Tgw(k)deltatw(k))
1481
1482	! GW formule 2
1483
1484	IF (dtimesub*tgw(i,k)<1.E-10) THEN
1485	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
1486	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1487	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
1488	tgw(i,k)*deltatw(i,k) )
1489	ELSE
1490	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i,k)))* &
1491	(ff(i)+dtke(i,k) - &
1492	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1493	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - &
1494	tgw(i,k)*deltatw(i,k) )
1495	END IF
1496
1497	dth(i, k) = deltatw(i, k)/ppi(i, k)
1498
1499	gg(i) = d_deltaqw(i, k)/dtimesub
1500
1501	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
1502	entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
1503	(death_rate(i)sigmaw(i)+detr(i,k))deltaqw(i,k)/(1.-sigmaw(i)))
1504	! cc
1505
1506	! cc nrlmd
1507	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1508	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1509	! cc
1510	END IF
1511	END DO
1512	END DO
1513
1514
1515	! Scale tendencies so that water vapour remains positive in w and x.
1516
1517	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
1518	d_deltaqw, sigmaw, d_sigmaw, alpha)
1519	!
1520	! Alpha_tot = Product of all the alpha's
1521	DO i = 1, klon
1522	IF (wk_adv(i)) THEN
1523	alpha_tot(i) = alpha_tot(i)*alpha(i)
1524	END IF
1525	END DO
1526
1527	! cc nrlmd
1528	! c print*,'alpha'
1529	! c do i=1,klon
1530	! c print*,alpha(i)
1531	! c end do
1532	! cc
1533	DO k = 1, klev
1534	DO i = 1, klon
1535	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1536	d_tenv(i, k) = alpha(i)*d_tenv(i, k)
1537	d_qe(i, k) = alpha(i)*d_qe(i, k)
1538	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1539	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1540	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1541	END IF
1542	END DO
1543	END DO
1544	DO i = 1, klon
1545	IF (wk_adv(i)) THEN
1546	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1547	END IF
1548	END DO
1549
1550	! Update large scale variables and wake variables
1551	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1552	! IM 060208 DO k = 1,kupper(i)
1553	DO k = 1, klev
1554	DO i = 1, klon
1555	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1556	dtls(i, k) = dtls(i, k) + d_tenv(i, k)
1557	dqls(i, k) = dqls(i, k) + d_qe(i, k)
1558	! cc nrlmd
1559	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1560	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1561	! cc
1562	END IF
1563	END DO
1564	END DO
1565	DO k = 1, klev
1566	DO i = 1, klon
1567	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1568	tenv(i, k) = tenv0(i, k) + dtls(i, k)
1569	qe(i, k) = qe0(i, k) + dqls(i, k)
1570	the(i, k) = tenv(i, k)/ppi(i, k)
1571	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1572	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1573	dth(i, k) = deltatw(i, k)/ppi(i, k)
1574	! c print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
1575	! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1576	END IF
1577	END DO
1578	END DO
1579	!
1580	DO i = 1, klon
1581	IF (wk_adv(i)) THEN
1582	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1583	d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
1584	END IF
1585	END DO
1586	!jyg<
1587	IF (iflag_wk_pop_dyn >= 1) THEN
1588	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1589	! Cumulatives
1590	DO i = 1, klon
1591	IF (wk_adv(i)) THEN
1592	d_sig_gen2(i) = d_sig_gen2(i) + d_sig_gen(i)
1593	d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
1594	d_sig_col2(i) = d_sig_col2(i) + d_sig_col(i)
1595	d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
1596	d_sig_bnd2(i) = d_sig_bnd2(i) + d_sig_bnd(i)
1597	END IF
1598	END DO
1599	! Bounds
1600	DO i = 1, klon
1601	IF (wk_adv(i)) THEN
1602	sigmaw_targ = max(sigmaw(i),sigmad)
1603	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
1604	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1605	sigmaw(i) = sigmaw_targ
1606	END IF
1607	END DO
1608	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1609	! Cumulatives
1610	DO i = 1, klon
1611	IF (wk_adv(i)) THEN
1612	wdens(i) = wdens(i) + d_wdens(i)
1613	d_wdens2(i) = d_wdens2(i) + d_wdens(i)
1614	d_dens_gen2(i) = d_dens_gen2(i) + d_dens_gen(i)
1615	d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
1616	d_dens_col2(i) = d_dens_col2(i) + d_dens_col(i)
1617	d_dens_bnd2(i) = d_dens_bnd2(i) + d_dens_bnd(i)
1618	END IF
1619	END DO
1620	! Bounds
1621	DO i = 1, klon
1622	IF (wk_adv(i)) THEN
1623	wdens_targ = max(wdens(i),wdensmin)
1624	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
1625	d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
1626	wdens(i) = wdens_targ
1627	END IF
1628	END DO
1629	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1630	! Cumulatives
1631	DO i = 1, klon
1632	IF (wk_adv(i)) THEN
1633	awdens(i) = awdens(i) + d_awdens(i)
1634	d_awdens2(i) = d_awdens2(i) + d_awdens(i)
1635	END IF
1636	END DO
1637	! Bounds
1638	DO i = 1, klon
1639	IF (wk_adv(i)) THEN
1640	wdens_targ = min( max(awdens(i),0.), wdens(i) )
1641	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
1642	d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
1643	awdens(i) = wdens_targ
1644	END IF
1645	END DO
1646	!
1647	IF (iflag_wk_pop_dyn >= 2) THEN
1648	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!!
1649	! Cumulatives
1650	DO i = 1, klon
1651	IF (wk_adv(i)) THEN
1652	d_adens_death2(i) = d_adens_death2(i) + d_adens_death(i)
1653	d_adens_icol2(i) = d_adens_icol2(i) + d_adens_icol(i)
1654	d_adens_acol2(i) = d_adens_acol2(i) + d_adens_acol(i)
1655	d_adens_bnd2(i) = d_adens_bnd2(i) + d_adens_bnd(i)
1656	END IF
1657	END DO
1658	! Bounds
1659	DO i = 1, klon
1660	IF (wk_adv(i)) THEN
1661	wdens_targ = min( max(awdens(i),0.), wdens(i) )
1662	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
1663	awdens(i) = wdens_targ
1664	END IF
1665	END DO
1666	!
1667	IF (iflag_wk_pop_dyn == 3) THEN
1668	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!!
1669	! Cumulatives
1670	DO i = 1, klon
1671	IF (wk_adv(i)) THEN
1672	asigmaw(i) = asigmaw(i) + d_asigmaw(i)
1673	d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i)
1674	d_asig_death2(i) = d_asig_death2(i) + d_asig_death(i)
1675	d_asig_spread2(i) = d_asig_spread2(i) + d_asig_spread(i)
1676	d_asig_iicol2(i) = d_asig_iicol2(i) + d_asig_iicol(i)
1677	d_asig_aicol2(i) = d_asig_aicol2(i) + d_asig_aicol(i)
1678	d_asig_bnd2(i) = d_asig_bnd2(i) + d_asig_bnd(i)
1679	END IF
1680	END DO
1681	! Bounds
1682	DO i = 1, klon
1683	IF (wk_adv(i)) THEN
1684	! asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad.
1685	!! sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
1686	sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i))
1687	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
1688	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
1689	asigmaw(i) = sigmaw_targ
1690	END IF
1691	END DO
1692
1693	#ifdef IOPHYS_WK
1694	IF (phys_sub) THEN
1695	CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens)
1696	CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens)
1697	CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw)
1698	CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw)
1699	!
1700	call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
1701	call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
1702	call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
1703	call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
1704	call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
1705	!
1706	call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
1707	call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
1708	call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
1709	call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
1710	call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
1711	!
1712	CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
1713	CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
1714	CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
1715	CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
1716	CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
1717	CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
1718	!
1719	CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
1720	CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
1721	CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
1722	CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
1723	CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
1724	CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
1725	ENDIF
1726	#endif
1727	ENDIF ! (iflag_wk_pop_dyn == 3)
1728	ENDIF ! (iflag_wk_pop_dyn >= 2)
1729	ENDIF ! (iflag_wk_pop_dyn >= 1)
1730
1731
1732
1733	Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
1734	dth, hw, rho, &
1735	ktop, wk_adv)
1736
1737
1738	! 5/ Set deltatw & deltaqw to 0 above kupper
1739
1740	DO k = 1, klev
1741	DO i = 1, klon
1742	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1743	deltatw(i, k) = 0.
1744	deltaqw(i, k) = 0.
1745	d_deltatw2(i,k) = -deltatw0(i,k)
1746	d_deltaqw2(i,k) = -deltaqw0(i,k)
1747	END IF
1748	END DO
1749	END DO
1750
1751
1752	! -------------Cstar computation---------------------------------
1753	DO i = 1, klon
1754	IF (wk_adv(i)) THEN !!! nrlmd
1755	sum_thu(i) = 0.
1756	sum_tu(i) = 0.
1757	sum_qu(i) = 0.
1758	sum_thvu(i) = 0.
1759	sum_dth(i) = 0.
1760	sum_dq(i) = 0.
1761	sum_rho(i) = 0.
1762	sum_dtdwn(i) = 0.
1763	sum_dqdwn(i) = 0.
1764
1765	av_thu(i) = 0.
1766	av_tu(i) = 0.
1767	av_qu(i) = 0.
1768	av_thvu(i) = 0.
1769	av_dth(i) = 0.
1770	av_dq(i) = 0.
1771	av_rho(i) = 0.
1772	av_dtdwn(i) = 0.
1773	av_dqdwn(i) = 0.
1774	END IF
1775	END DO
1776
1777	! Integrals (and wake top level number)
1778	! --------------------------------------
1779
1780	! Initialize sum_thvu to 1st level virt. pot. temp.
1781
1782	DO i = 1, klon
1783	IF (wk_adv(i)) THEN !!! nrlmd
1784	z(i) = 1.
1785	dz(i) = 1.
1786	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
1787	sum_dth(i) = 0.
1788	END IF
1789	END DO
1790
1791	DO k = 1, klev
1792	DO i = 1, klon
1793	IF (wk_adv(i)) THEN !!! nrlmd
1794	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
1795	IF (dz(i)>0) THEN
1796	z(i) = z(i) + dz(i)
1797	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1798	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1799	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1800	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
1801	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1802	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1803	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1804	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1805	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1806	END IF
1807	END IF
1808	END DO
1809	END DO
1810
1811	DO i = 1, klon
1812	IF (wk_adv(i)) THEN !!! nrlmd
1813	hw0(i) = z(i)
1814	END IF
1815	END DO
1816
1817
1818	! - WAPE and mean forcing computation
1819	! ---------------------------------------
1820
1821	! ---------------------------------------
1822
1823	! Means
1824
1825	DO i = 1, klon
1826	IF (wk_adv(i)) THEN !!! nrlmd
1827	av_thu(i) = sum_thu(i)/hw0(i)
1828	av_tu(i) = sum_tu(i)/hw0(i)
1829	av_qu(i) = sum_qu(i)/hw0(i)
1830	av_thvu(i) = sum_thvu(i)/hw0(i)
1831	av_dth(i) = sum_dth(i)/hw0(i)
1832	av_dq(i) = sum_dq(i)/hw0(i)
1833	av_rho(i) = sum_rho(i)/hw0(i)
1834	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1835	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1836
1837	wape(i) = -RGhw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
1838	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
1839	END IF
1840	END DO
1841
1842
1843	! Filter out bad wakes
1844
1845	DO k = 1, klev
1846	DO i = 1, klon
1847	IF (wk_adv(i)) THEN !!! nrlmd
1848	IF (wape(i)<0.) THEN
1849	deltatw(i, k) = 0.
1850	deltaqw(i, k) = 0.
1851	dth(i, k) = 0.
1852	d_deltatw2(i,k) = -deltatw0(i,k)
1853	d_deltaqw2(i,k) = -deltaqw0(i,k)
1854	END IF
1855	END IF
1856	END DO
1857	END DO
1858
1859	DO i = 1, klon
1860	IF (wk_adv(i)) THEN !!! nrlmd
1861	IF (wape(i)<0.) THEN
1862	wape(i) = 0.
1863	cstar(i) = 0.
1864	hw(i) = hwmin
1865	!jyg<
1866	!! sigmaw(i) = max(sigmad, sigd_con(i))
1867	sigmaw_targ = max(sigmad, sigd_con(i))
1868	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
1869	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1870	sigmaw(i) = sigmaw_targ
1871	!
1872	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
1873	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
1874	asigmaw(i) = sigmaw_targ
1875	!>jyg
1876	fip(i) = 0.
1877	gwake(i) = .FALSE.
1878	ELSE
1879	cstar(i) = starksqrt(2.wape(i))
1880	gwake(i) = .TRUE.
1881	END IF
1882	END IF
1883	END DO
1884	!
1885	! ------------------------------------------------------------------------
1886	!
1887	END DO ! isubstep end sub-timestep loop
1888
1889	!
1890	! ------------------------------------------------------------------------
1891	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1892	! ------------------------------------------------------------------------
1893	!
1894
1895	#ifdef IOPHYS_WK
1896	IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape)
1897	#endif
1898	IF (prt_level>=10) THEN
1899	PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
1900	sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
1901	ENDIF
1902
1903
1904	! ----------------------------------------------------------
1905	! Determine wake final state; recompute wape, cstar, ktop;
1906	! filter out bad wakes.
1907	! ----------------------------------------------------------
1908
1909	! 2.1 - Undisturbed area and Wake integrals
1910	! ---------------------------------------------------------
1911
1912	DO i = 1, klon
1913	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
1914	IF (ok_qx_qw(i)) THEN
1915	! cc
1916	z(i) = 0.
1917	sum_thu(i) = 0.
1918	sum_tu(i) = 0.
1919	sum_qu(i) = 0.
1920	sum_thvu(i) = 0.
1921	sum_dth(i) = 0.
1922	sum_half_dth(i) = 0.
1923	sum_dq(i) = 0.
1924	sum_rho(i) = 0.
1925	sum_dtdwn(i) = 0.
1926	sum_dqdwn(i) = 0.
1927
1928	av_thu(i) = 0.
1929	av_tu(i) = 0.
1930	av_qu(i) = 0.
1931	av_thvu(i) = 0.
1932	av_dth(i) = 0.
1933	av_dq(i) = 0.
1934	av_rho(i) = 0.
1935	av_dtdwn(i) = 0.
1936	av_dqdwn(i) = 0.
1937
1938	dthmin(i) = -delta_t_min
1939	END IF
1940	END DO
1941	! Potential temperatures and humidity
1942	! ----------------------------------------------------------
1943
1944	DO k = 1, klev
1945	DO i = 1, klon
1946	! cc nrlmd IF ( wk_adv(i)) THEN
1947	IF (ok_qx_qw(i)) THEN
1948	! cc
1949	rho(i, k) = p(i, k)/(RD*tenv(i,k))
1950	IF (k==1) THEN
1951	rhoh(i, k) = ph(i, k)/(RD*tenv(i,k))
1952	zhh(i, k) = 0
1953	ELSE
1954	rhoh(i, k) = ph(i, k)2./(RD(tenv(i,k)+tenv(i,k-1)))
1955	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
1956	END IF
1957	the(i, k) = tenv(i, k)/ppi(i, k)
1958	thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1959	tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i)
1960	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
1961	rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k)))
1962	dth(i, k) = deltatw(i, k)/ppi(i, k)
1963	END IF
1964	END DO
1965	END DO
1966
1967	! Integrals (and wake top level number)
1968	! -----------------------------------------------------------
1969
1970	! Initialize sum_thvu to 1st level virt. pot. temp.
1971
1972	DO i = 1, klon
1973	! cc nrlmd IF ( wk_adv(i)) THEN
1974	IF (ok_qx_qw(i)) THEN
1975	! cc
1976	z(i) = 1.
1977	dz(i) = 1.
1978	dz_half(i) = 1.
1979	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
1980	sum_dth(i) = 0.
1981	END IF
1982	END DO
1983
1984	DO k = 1, klev
1985	DO i = 1, klon
1986	! cc nrlmd IF ( wk_adv(i)) THEN
1987	IF (ok_qx_qw(i)) THEN
1988	! cc
1989	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
1990	dz_half(i) = -(amax1(ph(i,k+1),0.5(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)RG)
1991	IF (dz(i)>0) THEN
1992	z(i) = z(i) + dz(i)
1993	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1994	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1995	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1996	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
1997	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1998	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1999	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
2000	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
2001	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
2002	!
2003	dthmin(i) = min(dthmin(i), dth(i,k))
2004	END IF
2005	IF (dz_half(i)>0) THEN
2006	sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
2007	END IF
2008	END IF
2009	END DO
2010	END DO
2011
2012	DO i = 1, klon
2013	! cc nrlmd IF ( wk_adv(i)) THEN
2014	IF (ok_qx_qw(i)) THEN
2015	! cc
2016	hw0(i) = z(i)
2017	END IF
2018	END DO
2019
2020	! - WAPE and mean forcing computation
2021	! -------------------------------------------------------------
2022
2023	! Means
2024
2025	DO i = 1, klon
2026	! cc nrlmd IF ( wk_adv(i)) THEN
2027	IF (ok_qx_qw(i)) THEN
2028	! cc
2029	av_thu(i) = sum_thu(i)/hw0(i)
2030	av_tu(i) = sum_tu(i)/hw0(i)
2031	av_qu(i) = sum_qu(i)/hw0(i)
2032	av_thvu(i) = sum_thvu(i)/hw0(i)
2033	av_dth(i) = sum_dth(i)/hw0(i)
2034	av_dq(i) = sum_dq(i)/hw0(i)
2035	av_rho(i) = sum_rho(i)/hw0(i)
2036	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
2037	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
2038
2039	wape2(i) = -RGhw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
2040	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
2041	END IF
2042	END DO
2043	#ifdef IOPHYS_WK
2044	IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2)
2045	#endif
2046
2047
2048	! Prognostic variable update
2049	! ------------------------------------------------------------
2050
2051	! Filter out bad wakes
2052
2053	IF (iflag_wk_check_trgl>=1) THEN
2054	! Check triangular shape of dth profile
2055	DO i = 1, klon
2056	IF (ok_qx_qw(i)) THEN
2057	!! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
2058	!! print ,'wake, 2.sum_dth(i)/(hw0(i)*dthmin(i)) ', &
2059	!! 2.sum_dth(i)/(hw0(i)dthmin(i))
2060	!! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
2061	!! sum_half_dth(i), sum_dth(i)
2062	IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
2063	wape2(i) = -1.
2064	!! print *,'wake, rej 1'
2065	ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.sum_dth(i)/(hw0(i)dthmin(i)) - 1.) > 0.5) THEN
2066	wape2(i) = -1.
2067	!! print *,'wake, rej 2'
2068	ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
2069	wape2(i) = -1.
2070	!! print *,'wake, rej 3'
2071	END IF
2072	END IF
2073	END DO
2074	END IF
2075	#ifdef IOPHYS_WK
2076	IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)
2077	#endif
2078
2079
2080	DO k = 1, klev
2081	DO i = 1, klon
2082	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
2083	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
2084	! cc
2085	deltatw(i, k) = 0.
2086	deltaqw(i, k) = 0.
2087	dth(i, k) = 0.
2088	d_deltatw2(i,k) = -deltatw0(i,k)
2089	d_deltaqw2(i,k) = -deltaqw0(i,k)
2090	END IF
2091	END DO
2092	END DO
2093
2094
2095	DO i = 1, klon
2096	! cc nrlmd IF ( wk_adv(i)) THEN
2097	IF (ok_qx_qw(i)) THEN
2098	! cc
2099	IF (wape2(i)<0.) THEN
2100	wape2(i) = 0.
2101	cstar2(i) = 0.
2102	hw(i) = hwmin
2103	!jyg<
2104	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
2105	sigmaw_targ = max(sigmad, sigd_con(i))
2106	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
2107	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2108	sigmaw(i) = sigmaw_targ
2109	!
2110	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
2111	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
2112	asigmaw(i) = sigmaw_targ
2113	!>jyg
2114	fip(i) = 0.
2115	gwake(i) = .FALSE.
2116	ELSE
2117	IF (prt_level>=10) PRINT *, 'wape2>0'
2118	cstar2(i) = starksqrt(2.wape2(i))
2119	gwake(i) = .TRUE.
2120	END IF
2121	#ifdef IOPHYS_WK
2122	IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2)
2123	#endif
2124	END IF ! (ok_qx_qw(i))
2125	END DO
2126
2127	DO i = 1, klon
2128	! cc nrlmd IF ( wk_adv(i)) THEN
2129	IF (ok_qx_qw(i)) THEN
2130	! cc
2131	ktopw(i) = ktop(i)
2132	END IF
2133	END DO
2134
2135	DO i = 1, klon
2136	! cc nrlmd IF ( wk_adv(i)) THEN
2137	IF (ok_qx_qw(i)) THEN
2138	! cc
2139	IF (ktopw(i)>0 .AND. gwake(i)) THEN
2140
2141	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
2142	! cc heff = 600.
2143	! Utilisation de la hauteur hw
2144	! c heff = 0.7*hw
2145	heff(i) = hw(i)
2146
2147	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
2148	sqrt(sigmaw(i)wdens(i)3.14)
2149	fip(i) = alpk*fip(i)
2150	! jyg2
2151	ELSE
2152	fip(i) = 0.
2153	END IF
2154	END IF
2155	END DO
2156	IF (iflag_wk_pop_dyn >= 3) THEN
2157	#ifdef IOPHYS_WK
2158	IF (.not.phys_sub) THEN
2159	CALL iophys_ecrit('fip',1,'fip','J/kg',fip)
2160	CALL iophys_ecrit('hw',1,'hw','J/kg',hw)
2161	CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop)
2162	CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens)
2163	CALL iophys_ecrit('awdens',1,'awdens','m',awdens)
2164	CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw)
2165	CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw)
2166	!
2167	CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk)
2168	CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk)
2169	CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk)
2170	!
2171	call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
2172	call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
2173	call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
2174	call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
2175	call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
2176	!
2177	call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
2178	call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
2179	call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
2180	call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
2181	call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
2182	!
2183	CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
2184	CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
2185	CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
2186	CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
2187	CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
2188	CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
2189	!
2190	CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
2191	CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
2192	CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
2193	CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
2194	CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
2195	CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
2196	ENDIF ! (.not.phys_sub)
2197	#endif
2198	ENDIF ! (iflag_wk_pop_dyn >= 3)
2199	! Limitation de sigmaw
2200
2201	! cc nrlmd
2202	! DO i=1,klon
2203	! IF (OK_qx_qw(i)) THEN
2204	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
2205	! ENDIF
2206	! ENDDO
2207	! cc
2208
2209	!jyg<
2210	IF (iflag_wk_pop_dyn >= 1) THEN
2211	DO i = 1, klon
2212	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2213	.NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold)
2214	!! .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
2215	ENDDO
2216	ELSE ! (iflag_wk_pop_dyn >= 1)
2217	DO i = 1, klon
2218	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2219	.NOT. ok_qx_qw(i)
2220	ENDDO
2221	ENDIF ! (iflag_wk_pop_dyn >= 1)
2222	!>jyg
2223
2224	DO k = 1, klev
2225	DO i = 1, klon
2226	!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2227	!!jyg .NOT. ok_qx_qw(i)) THEN
2228	IF (kill_wake(i)) THEN
2229	! cc
2230	dtls(i, k) = 0.
2231	dqls(i, k) = 0.
2232	deltatw(i, k) = 0.
2233	deltaqw(i, k) = 0.
2234	d_deltatw2(i,k) = -deltatw0(i,k)
2235	d_deltaqw2(i,k) = -deltaqw0(i,k)
2236	END IF ! (kill_wake(i))
2237	END DO
2238	END DO
2239
2240	DO i = 1, klon
2241	!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2242	!!jyg .NOT. ok_qx_qw(i)) THEN
2243	IF (kill_wake(i)) THEN
2244	ktopw(i) = 0
2245	wape(i) = 0.
2246	cstar(i) = 0.
2247	!!jyg Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes
2248	!! hw(i) = hwmin !jyg
2249	!! sigmaw(i) = sigmad !jyg
2250	hw(i) = 0. !jyg
2251	fip(i) = 0.
2252	!
2253	!! sigmaw(i) = 0. !jyg
2254	sigmaw_targ = 0.
2255	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
2256	!! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2257	d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i) ! _in = correction jyg 20220124
2258	sigmaw(i) = sigmaw_targ
2259	!
2260	IF (iflag_wk_pop_dyn >= 3) THEN
2261	sigmaw_targ = 0.
2262	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
2263	!! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2264	d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i) ! _in = correction jyg 20220124
2265	asigmaw(i) = sigmaw_targ
2266	ELSE
2267	asigmaw(i) = 0.
2268	ENDIF ! (iflag_wk_pop_dyn >= 3)
2269	!
2270	IF (iflag_wk_pop_dyn >= 1) THEN
2271	!! awdens(i) = 0.
2272	!! wdens(i) = 0.
2273	wdens_targ = 0.
2274	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
2275	!! d_wdens2(i) = wdens_targ - wdens(i)
2276	d_wdens2(i) = wdens_targ - wdens_in(i) ! jyg 20220916
2277	wdens(i) = wdens_targ
2278	wdens_targ = 0.
2279	!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
2280	IF (iflag_wk_pop_dyn >= 2) THEN
2281	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
2282	ENDIF ! (iflag_wk_pop_dyn >= 2)
2283	!! d_awdens2(i) = wdens_targ - awdens(i)
2284	d_awdens2(i) = wdens_targ - awdens_in(i) ! jyg 20220916
2285	awdens(i) = wdens_targ
2286	!! IF (iflag_wk_pop_dyn == 2) THEN
2287	!! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
2288	!! ENDIF ! (iflag_wk_pop_dyn == 2)
2289	ENDIF ! (iflag_wk_pop_dyn >= 1)
2290	ELSE ! (kill_wake(i))
2291	wape(i) = wape2(i)
2292	cstar(i) = cstar2(i)
2293	END IF ! (kill_wake(i))
2294	! c print*,'wape wape2 ktopw OK_qx_qw =',
2295	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
2296	END DO
2297
2298	IF (prt_level>=10) THEN
2299	PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
2300	wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
2301	ENDIF
2302	#ifdef IOPHYS_WK
2303	IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)
2304	#endif
2305
2306
2307	! -----------------------------------------------------------------
2308	! Get back to tendencies per second
2309
2310	DO k = 1, klev
2311	DO i = 1, klon
2312
2313	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
2314	!jyg<
2315	!! IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
2316	IF (ok_qx_qw(i)) THEN
2317	!>jyg
2318	! cc
2319	dtls(i, k) = dtls(i, k)/dtime
2320	dqls(i, k) = dqls(i, k)/dtime
2321	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
2322	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
2323	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
2324	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
2325	! c $ ,death_rate(i)*sigmaw(i)
2326	END IF
2327	END DO
2328	END DO
2329	!jyg<
2330	IF (iflag_wk_pop_dyn >= 1) THEN
2331	DO i = 1, klon
2332	IF (ok_qx_qw(i)) THEN
2333	d_sig_gen2(i) = d_sig_gen2(i)/dtime
2334	d_sig_death2(i) = d_sig_death2(i)/dtime
2335	d_sig_col2(i) = d_sig_col2(i)/dtime
2336	d_sig_spread2(i) = d_sig_spread2(i)/dtime
2337	d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
2338	d_sigmaw2(i) = d_sigmaw2(i)/dtime
2339	!
2340	d_dens_gen2(i) = d_dens_gen2(i)/dtime
2341	d_dens_death2(i) = d_dens_death2(i)/dtime
2342	d_dens_col2(i) = d_dens_col2(i)/dtime
2343	d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
2344	d_awdens2(i) = d_awdens2(i)/dtime
2345	d_wdens2(i) = d_wdens2(i)/dtime
2346	ENDIF
2347	ENDDO
2348	IF (iflag_wk_pop_dyn >= 2) THEN
2349	DO i = 1, klon
2350	IF (ok_qx_qw(i)) THEN
2351	d_adens_death2(i) = d_adens_death2(i)/dtime
2352	d_adens_icol2(i) = d_adens_icol2(i)/dtime
2353	d_adens_acol2(i) = d_adens_acol2(i)/dtime
2354	d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
2355	ENDIF
2356	ENDDO
2357	IF (iflag_wk_pop_dyn == 3) THEN
2358	DO i = 1, klon
2359	IF (ok_qx_qw(i)) THEN
2360	d_asig_death2(i) = d_asig_death2(i)/dtime
2361	d_asig_iicol2(i) = d_asig_iicol2(i)/dtime
2362	d_asig_aicol2(i) = d_asig_aicol2(i)/dtime
2363	d_asig_spread2(i) = d_asig_spread2(i)/dtime
2364	d_asig_bnd2(i) = d_asig_bnd2(i)/dtime
2365	ENDIF
2366	ENDDO
2367	ENDIF ! (iflag_wk_pop_dyn == 3)
2368	ENDIF ! (iflag_wk_pop_dyn >= 2)
2369	ENDIF ! (iflag_wk_pop_dyn >= 1)
2370
2371	!>jyg
2372
2373	RETURN
2374	END SUBROUTINE wake
2375
2376	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
2377	d_deltaqw, sigmaw, d_sigmaw, alpha)
2378	! ------------------------------------------------------
2379	! Dtermination du coefficient alpha tel que les tendances
2380	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
2381	! a une humidite positive dans la zone (x) et dans la zone (w).
2382	! ------------------------------------------------------
2383	IMPLICIT NONE
2384
2385	! Input
2386	REAL qe(nlon, nl), d_qe(nlon, nl)
2387	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
2388	REAL sigmaw(nlon), d_sigmaw(nlon)
2389	LOGICAL wk_adv(nlon)
2390	INTEGER nl, nlon
2391	! Output
2392	REAL alpha(nlon)
2393	! Internal variables
2394	REAL zeta(nlon, nl)
2395	REAL alpha1(nlon)
2396	REAL x, a, b, c, discrim
2397	REAL epsilon_loc
2398	INTEGER i,k
2399
2400	DO k = 1, nl
2401	DO i = 1, nlon
2402	IF (wk_adv(i)) THEN
2403	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
2404	zeta(i, k) = 0.
2405	ELSE
2406	zeta(i, k) = 1.
2407	END IF
2408	END IF
2409	END DO
2410	DO i = 1, nlon
2411	IF (wk_adv(i)) THEN
2412	x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
2413	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i) &
2414	(deltaqw(i,k)+d_deltaqw(i,k))
2415	a = -d_sigmaw(i)*d_deltaqw(i, k)
2416	b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
2417	deltaqw(i, k)*d_sigmaw(i)
2418	c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc
2419	discrim = bb - 4.a*c
2420	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
2421	IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap
2422	alpha1(i) = 1.
2423	ELSE
2424	IF (x>=0.) THEN
2425	alpha1(i) = 1.
2426	ELSE
2427	IF (a>0.) THEN
2428	alpha1(i) = 0.9min( (2.c)/(-b+sqrt(discrim)), &
2429	(-b+sqrt(discrim))/(2.*a) )
2430	ELSE IF (a==0.) THEN
2431	alpha1(i) = 0.9*(-c/b)
2432	ELSE
2433	! print*,'a,b,c discrim',a,b,c discrim
2434	alpha1(i) = 0.9max( (2.c)/(-b+sqrt(discrim)), &
2435	(-b+sqrt(discrim))/(2.*a))
2436	END IF
2437	END IF
2438	END IF
2439	alpha(i) = min(alpha(i), alpha1(i))
2440	END IF
2441	END DO
2442	END DO
2443
2444	RETURN
2445	END SUBROUTINE wake_vec_modulation
2446
2447
2448
2449	SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
2450	dth, hw_, rho, &
2451	ktop, wk_adv)
2452
2453	USE lmdz_wake_ini , ONLY : wk_pupper
2454	USE lmdz_wake_ini , ONLY : RG
2455	USE lmdz_wake_ini , ONLY : hwmin
2456	IMPLICIT NONE
2457
2458	INTEGER, INTENT(IN) :: klon,klev
2459	REAL, INTENT(IN), DIMENSION (klon,klev+1) :: ph, p
2460	REAL, INTENT(IN), DIMENSION (klon,klev+1) :: rho
2461	LOGICAL, INTENT(IN), DIMENSION (klon) :: wk_adv
2462	REAL, INTENT(IN), DIMENSION (klon,klev+1) :: dth
2463
2464	REAL, INTENT(OUT), DIMENSION (klon) :: hw_
2465	REAL, INTENT(OUT), DIMENSION (klon) :: ptop
2466	INTEGER, INTENT(OUT), DIMENSION (klon) :: Ktop
2467	REAL, INTENT(OUT), DIMENSION (klon) :: pupper
2468	INTEGER, INTENT(OUT), DIMENSION (klon) :: kupper
2469	INTEGER :: i,k
2470	REAL :: delta_t_min
2471
2472
2473	REAL, DIMENSION (klon) :: dthmin
2474	REAL, DIMENSION (klon) :: ptop_provis,ptop_new
2475	REAL, DIMENSION (klon) :: z, dz
2476	REAL, DIMENSION (klon) :: sum_dth
2477
2478	!INTEGER, SAVE :: compte=0
2479
2480	! LJYF : a priori z, dz sum_dth sont aussi des variables internes
2481	! Les eliminer apres verification convergence numerique
2482
2483	!compte=compte+1
2484	!print*,'compte=',compte
2485
2486	delta_t_min = 0.2
2487
2488	! Determine Ptop from buoyancy integral
2489	! ---------------------------------------
2490
2491	! - 1/ Pressure of the level where dth changes sign.
2492	!print*,'WAKE LJYF'
2493
2494	DO i = 1, klon
2495	ptop_provis(i) = ph(i, 1)
2496	END DO
2497
2498	DO k = 2, klev
2499	DO i = 1, klon
2500	! if (compte==92) then
2501	! print*,'debug xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
2502	! print*,'debug i =',i
2503	! print*,'debug k =',k
2504	! print*,'debug wk_adv(i) =',wk_adv(i)
2505	! print*,'debug ptop_provis(i) =',ptop_provis(i)
2506	! print*,'debug ph(i,1) =',ph(i,1)
2507	! print*,'debug dth(i,k) =',dth(i,k)
2508	! print*,'debug delta_t_min =',delta_t_min
2509	! print*,'debug p(i,k-1) =',p(i,k-1)
2510	! print*,'debug dth(i,k-1) =',dth(i,k-1)
2511	! print*,'debug p(i,k) =',p(i,k)
2512	! endif
2513	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
2514	! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2515	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2516	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
2517	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
2518	END IF
2519	END DO
2520	END DO
2521
2522	! - 2/ dth integral
2523
2524	DO i = 1, klon
2525	IF (wk_adv(i)) THEN !!! nrlmd
2526	sum_dth(i) = 0.
2527	dthmin(i) = -delta_t_min
2528	z(i) = 0.
2529	END IF
2530	END DO
2531
2532	DO k = 1, klev
2533	DO i = 1, klon
2534	IF (wk_adv(i)) THEN
2535	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
2536	IF (dz(i)>0) THEN
2537	z(i) = z(i) + dz(i)
2538	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
2539	dthmin(i) = amin1(dthmin(i), dth(i,k))
2540	END IF
2541	END IF
2542	END DO
2543	END DO
2544
2545	! - 3/ height of triangle with area= sum_dth and base = dthmin
2546
2547	DO i = 1, klon
2548	IF (wk_adv(i)) THEN
2549	hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
2550	hw_(i) = amax1(hwmin, hw_(i))
2551	END IF
2552	END DO
2553
2554	! - 4/ now, get Ptop
2555
2556	DO i = 1, klon
2557	IF (wk_adv(i)) THEN !!! nrlmd
2558	ktop(i) = 0
2559	z(i) = 0.
2560	END IF
2561	END DO
2562
2563	DO k = 1, klev
2564	DO i = 1, klon
2565	IF (wk_adv(i)) THEN
2566	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i))
2567	IF (dz(i)>0) THEN
2568	z(i) = z(i) + dz(i)
2569	ptop(i) = ph(i, k) - rho(i, k)RGdz(i)
2570	ktop(i) = k
2571	END IF
2572	END IF
2573	END DO
2574	END DO
2575
2576	! 4.5/Correct ktop and ptop
2577
2578	DO i = 1, klon
2579	ptop_new(i) = ptop(i)
2580	END DO
2581
2582	DO k = klev, 2, -1
2583	DO i = 1, klon
2584	! IM v3JYG; IF (k .GE. ktop(i)
2585	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
2586	! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2587	dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2588	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
2589	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
2590	END IF
2591	END DO
2592	END DO
2593
2594
2595	DO i = 1, klon
2596	ptop(i) = ptop_new(i)
2597	END DO
2598
2599	DO k = klev, 1, -1
2600	DO i = 1, klon
2601	IF (wk_adv(i)) THEN !!! nrlmd
2602	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
2603	END IF
2604	END DO
2605	END DO
2606
2607	! IF (prt_level>=10) THEN
2608	! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
2609	! ENDIF
2610
2611	kupper = 0
2612
2613	IF (wk_pupper<1.) THEN
2614	! Choose an integration bound well above wake top
2615	! -----------------------------------------------------------------
2616
2617	! Pupper = 50000. ! melting level
2618	! Pupper = 60000.
2619	! Pupper = 80000. ! essais pour case_e
2620	DO i = 1, klon
2621	! pupper(i) = 0.6*ph(i, 1)
2622	pupper(i) = wk_pupper*ph(i, 1)
2623	pupper(i) = max(pupper(i), 45000.)
2624	! cc Pupper(i) = 60000.
2625	END DO
2626
2627	ELSE
2628
2629	DO i=1, klon
2630	! pupper(i) = wk_pupperptop(i)+(1.-wk_pupper)ph(i, 1)
2631	pupper(i) = min( wk_pupperptop(i)+(1.-wk_pupper)ph(i, 1) , ptop(i)-5000.)
2632	END DO
2633	END IF
2634
2635	! -5/ Determination de kupper
2636
2637	DO k = klev, 1, -1
2638	DO i = 1, klon
2639	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
2640	END DO
2641	END DO
2642
2643	! On evite kupper = 1 et kupper = klev
2644	DO i = 1, klon
2645	kupper(i) = max(kupper(i), 2)
2646	kupper(i) = min(kupper(i), klev-1)
2647	END DO
2648	RETURN
2649	END SUBROUTINE pkupper
2650
2651
2652	SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
2653	wdensmin, &
2654	dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
2655	d_awdens, d_wdens, d_sigmaw, &
2656	iflag_wk_act, wk_adv, cin, wape, &
2657	drdt, &
2658	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
2659	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
2660	d_wdens_targ, d_sigmaw_targ)
2661
2662
2663	USE lmdz_wake_ini , ONLY : wake_ini
2664	USE lmdz_wake_ini , ONLY : prt_level,RG
2665	USE lmdz_wake_ini , ONLY : stark, wdens_ref
2666	USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
2667	!! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
2668	USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
2669	USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
2670
2671	IMPLICIT NONE
2672
2673	INTEGER, INTENT(IN) :: klon,klev
2674	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
2675	REAL, INTENT(IN) :: dtime
2676	REAL, INTENT(IN) :: dtimesub
2677	REAL, INTENT(IN) :: wdensmin
2678	REAL, DIMENSION (klon), INTENT(IN) :: wgen
2679	REAL, DIMENSION (klon), INTENT(IN) :: wdens
2680	REAL, DIMENSION (klon), INTENT(IN) :: awdens
2681	REAL, DIMENSION (klon), INTENT(IN) :: sigmaw
2682	REAL, DIMENSION (klon), INTENT(IN) :: cstar
2683	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape
2684	REAL, DIMENSION (klon), INTENT(IN) :: f_shear
2685	INTEGER, INTENT(IN) :: iflag_wk_act
2686
2687
2688	!
2689
2690	! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
2691	! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
2692	REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk
2693	REAL, DIMENSION (klon), INTENT(OUT) :: gfl
2694	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_awdens, d_wdens
2695	REAL, DIMENSION (klon), INTENT(OUT) :: drdt
2696	! Some components of the tendencies of state variables
2697	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
2698	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_spread
2699	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
2700	REAL, INTENT(OUT) :: d_wdens_targ, d_sigmaw_targ
2701
2702
2703	REAL :: delta_t_min
2704	INTEGER :: nsub
2705	INTEGER :: i, k
2706	REAL :: wdens0
2707	! IM 080208
2708	LOGICAL, DIMENSION (klon) :: gwake
2709
2710	! Variables liees a la dynamique de population
2711	REAL, DIMENSION(klon) :: act
2712	REAL, DIMENSION(klon) :: tau_wk_inv
2713	REAL, DIMENSION(klon) :: wape1_act, wape2_act
2714	LOGICAL, DIMENSION (klon) :: kill_wake
2715	REAL :: drdt_pos
2716	REAL :: tau_wk_inv_min
2717
2718
2719
2720	IF (iflag_wk_act == 0) THEN
2721	act(:) = 0.
2722	ELSEIF (iflag_wk_act == 1) THEN
2723	act(:) = 1.
2724	ELSEIF (iflag_wk_act ==2) THEN
2725	DO i = 1, klon
2726	IF (wk_adv(i)) THEN
2727	wape1_act(i) = abs(cin(i))
2728	wape2_act(i) = 2.*wape1_act(i) + 1.
2729	act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
2730	ENDIF ! (wk_adv(i))
2731	ENDDO
2732	ENDIF ! (iflag_wk_act ==2)
2733
2734	DO i = 1, klon
2735	IF (wk_adv(i)) THEN
2736	rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
2737	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
2738	END IF
2739	END DO
2740
2741	DO i = 1, klon
2742	IF (wk_adv(i)) THEN
2743	!! tau_wk(i) = max(rad_wk(i)/(3.cstar(i))((cstar(i)/cstart)**1.5 - 1), 100.)
2744	tau_wk_inv(i) = max( (3.cstar(i))/(rad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
2745	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
2746	drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
2747	(1 + 2f_shear(i)(2.sigmaw(i)-aa0wdens(i)) - 2.*sigmaw(i))
2748	!! (1 - 2sigmaw(i)(1.-f_shear(i)))
2749	drdt_pos=max(drdt(i),0.)
2750
2751	!! d_wdens(i) = ( wgen(i)(1.+2.(sigmaw(i)-sigmad)) &
2752	!! - wdens(i)*tau_wk_inv_min &
2753	!! - 2.gfl(i)wdens(i)Cstar(i) )dtimesub
2754	!jyg+mlt<
2755	d_awdens(i) = ( wgen(i) - (1./tau_cv)(awdens(i) - act(i)wdens(i)) )*dtimesub
2756	d_dens_gen(i) = wgen(i)
2757	d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
2758	d_dens_col(i) = -2.wdens(i)gfl(i)*drdt_pos
2759	d_dens_gen(i) = d_dens_gen(i)*dtimesub
2760	d_dens_death(i) = d_dens_death(i)*dtimesub
2761	d_dens_col(i) = d_dens_col(i)*dtimesub
2762
2763	d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
2764	!! d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min - &
2765	!! 2.wdens(i)gfl(i)drdt_pos )dtimesub
2766	!>jyg+mlt
2767	!
2768	!jyg<
2769	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
2770	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
2771	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
2772	d_wdens(i) = d_wdens_targ
2773	!! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
2774	!>jyg
2775
2776	!jyg+mlt<
2777	!! d_sigmaw(i) = ( (1.-2f_shear(i)sigmaw(i))(gfl(i)Cstar(i)+wgen(i)*sigmad/wdens(i)) &
2778	!! + 2.f_shear(i)wgen(i)sigmaw(i)*2/wdens(i) &
2779	!! - sigmaw(i)tau_wk_inv_min )dtimesub
2780	d_sig_gen(i) = wgen(i)*aa0
2781	d_sig_death(i) = - sigmaw(i)(1.-awdens(i)/wdens(i))tau_wk_inv_min
2782	!!
2783
2784	d_sig_col(i) = - 2f_shear(i)sigmaw(i)gfl(i)drdt_pos
2785	d_sig_col(i) = - 2f_shear(i)(2.sigmaw(i)-wdens(i)aa0)gfl(i)drdt_pos
2786	d_sig_spread(i) = gfl(i)*cstar(i)
2787	d_sig_gen(i) = d_sig_gen(i)*dtimesub
2788	d_sig_death(i) = d_sig_death(i)*dtimesub
2789	d_sig_col(i) = d_sig_col(i)*dtimesub
2790	d_sig_spread(i) = d_sig_spread(i)*dtimesub
2791	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
2792	!>jyg+mlt
2793	!
2794	!jyg<
2795	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
2796	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
2797	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
2798	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
2799	d_sigmaw(i) = d_sigmaw_targ
2800	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
2801	!>jyg
2802	ENDIF
2803	ENDDO
2804
2805	IF (prt_level >= 10) THEN
2806	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
2807	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
2808	print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
2809	wdens(1), awdens(1), act(1), d_awdens(1)
2810	print ,'wake, wgen, -(wdens-awdens)tau_wk_inv, -2.wdensgfl*drdt_pos, d_wdens ', &
2811	wgen(1), -(wdens(1)-awdens(1))tau_wk_inv(1), -2.wdens(1)gfl(1)drdt_pos, d_wdens(1)
2812	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
2813	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
2814	ENDIF
2815
2816	RETURN
2817	END SUBROUTINE wake_popdyn_1
2818
2819	SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
2820	wdensmin, &
2821	sigmaw, wdens, awdens, & !! states variables
2822	gfl, cstar, cin, wape, rad_wk, &
2823	d_sigmaw, d_wdens, d_awdens, & !! tendences
2824	cont_fact, &
2825	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
2826	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
2827	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
2828
2829
2830
2831	USE lmdz_wake_ini , ONLY : wake_ini
2832	USE lmdz_wake_ini , ONLY : prt_level,RG
2833	USE lmdz_wake_ini , ONLY : stark, wdens_ref
2834	USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
2835	!! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
2836	USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
2837	USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
2838
2839	IMPLICIT NONE
2840
2841	INTEGER, INTENT(IN) :: klon,klev
2842	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
2843	REAL, INTENT(IN) :: dtimesub
2844	REAL, INTENT(IN) :: wdensmin
2845	REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes
2846	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes
2847	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area
2848	REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area
2849	REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes
2850	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre
2851
2852	!
2853	REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = wake radius
2854	REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front lenght per unit area
2855	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens
2856	REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact !! RM facteur de contact = 2 pi * rad * C*
2857	! Some components of the tendencies of state variables
2858	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
2859	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
2860	REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
2861
2862
2863	!! internal variables
2864
2865	INTEGER :: i, k
2866	REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes
2867	REAL :: tau_wk_inv_min
2868	REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes
2869	REAL :: d_wdens_targ, d_sigmaw_targ
2870
2871
2872	!! Equations
2873	!! dD/dt = B - (D-A)/tau - f D^2
2874	!! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2
2875	!! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0)
2876	!!
2877	!! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma))
2878
2879
2880	DO i = 1, klon
2881	IF (wk_adv(i)) THEN
2882	rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
2883	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
2884	END IF
2885	END DO
2886
2887
2888	DO i = 1, klon
2889	IF (wk_adv(i)) THEN
2890	!! tau_wk(i) = max(rad_wk(i)/(3.cstar(i))((cstar(i)/cstart)**1.5 - 1), 100.)
2891	tau_wk_inv(i) = max( (3.cstar(i))/(rad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
2892	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
2893	tau_prime(i) = tau_cv
2894	!! cont_fact(i) = 2.(wgen(i)(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / &
2895	!! (2.(wdens(i)-awdens(i))2(2.sigmaw(i)/wdens(i) - aa0) + wdens(i)(1.-2.*sigmaw(i)))
2896	!! cont_fact(i) = 2.3.14rad_wk(i)*cstar(i) ! bug
2897	!! cont_fact(i) = 4.3.14rad_wk(i)*cstar(i)
2898	cont_fact(i) = 2.gfl(i)cstar(i)/wdens(i)
2899
2900	d_sig_gen(i) = wgen(i)*aa0
2901	d_sig_death(i) = - sigmaw(i)(1.-awdens(i)/wdens(i))tau_wk_inv_min
2902	d_sig_col(i) = - cont_fact(i)(wdens(i)-awdens(i))2(2.*sigmaw(i)/wdens(i)-aa0)
2903	d_sig_spread(i) = gfl(i)*cstar(i)
2904	!
2905	d_sig_gen(i) = d_sig_gen(i)*dtimesub
2906	d_sig_death(i) = d_sig_death(i)*dtimesub
2907	d_sig_col(i) = d_sig_col(i)*dtimesub
2908	d_sig_spread(i) = d_sig_spread(i)*dtimesub
2909	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
2910
2911
2912	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
2913	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
2914	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
2915	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
2916	d_sigmaw(i) = d_sigmaw_targ
2917	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
2918
2919
2920	d_dens_gen(i) = wgen(i)
2921	d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
2922	d_dens_col(i) = - cont_fact(i)wdens(i)*2
2923	!
2924	d_dens_gen(i) = d_dens_gen(i)*dtimesub
2925	d_dens_death(i) = d_dens_death(i)*dtimesub
2926	d_dens_col(i) = d_dens_col(i)*dtimesub
2927	d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
2928
2929
2930	d_adens_death(i) = -awdens(i)/tau_prime(i)
2931	d_adens_icol(i) = cont_fact(i)(wdens(i)-awdens(i))*2
2932	d_adens_acol(i) = - cont_fact(i)awdens(i)*2
2933	!
2934	d_adens_death(i) = d_adens_death(i)*dtimesub
2935	d_adens_icol(i) = d_adens_icol(i)*dtimesub
2936	d_adens_acol(i) = d_adens_acol(i)*dtimesub
2937	d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)
2938
2939	!!
2940	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
2941	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
2942	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
2943	d_wdens(i) = d_wdens_targ
2944
2945	d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
2946	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
2947	d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
2948	d_awdens(i) = d_wdens_targ
2949
2950
2951
2952	ENDIF
2953	ENDDO
2954
2955	IF (prt_level >= 10) THEN
2956	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', &
2957	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1)
2958	print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
2959	wdens(1), awdens(1), d_awdens(1)
2960	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
2961	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
2962	ENDIF
2963	sigmaw=sigmaw+d_sigmaw
2964	wdens=wdens+d_wdens
2965	awdens=awdens+d_awdens
2966
2967	RETURN
2968	END SUBROUTINE wake_popdyn_2
2969
2970	SUBROUTINE wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
2971	wdensmin, &
2972	sigmaw, asigmaw, wdens, awdens, & !! state variables
2973	gfl, agfl, cstar, cin, wape, &
2974	rad_wk, arad_wk, irad_wk, &
2975	d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies
2976	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
2977	d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
2978	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
2979	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
2980
2981
2982
2983	USE lmdz_wake_ini , ONLY : wake_ini
2984	USE lmdz_wake_ini , ONLY : prt_level,RG
2985	USE lmdz_wake_ini , ONLY : stark, wdens_ref
2986	USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
2987	!! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
2988	USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
2989	USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
2990	USE lmdz_wake_ini , ONLY : smallestreal
2991
2992	IMPLICIT NONE
2993
2994	INTEGER, INTENT(IN) :: klon,klev
2995	LOGICAL, INTENT(IN) :: phys_sub
2996	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
2997	REAL, INTENT(IN) :: dtimesub
2998	REAL, INTENT(IN) :: wdensmin
2999	REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes
3000	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes
3001	REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw !! sigma = fractional area of active wakes
3002	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area
3003	REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area
3004	REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes
3005	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre
3006
3007	!
3008	REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = mean wake radius
3009	REAL, DIMENSION (klon), INTENT(OUT) :: arad_wk !! r_A = wake radius of active wakes
3010	REAL, DIMENSION (klon), INTENT(OUT) :: irad_wk !! r_I = wake radius of inactive wakes
3011	REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front length per unit area
3012	REAL, DIMENSION (klon), INTENT(OUT) :: agfl !! LgA = gust front length of active wakes
3013	!! per unit area
3014	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_asigmaw, d_wdens, d_awdens
3015	! Some components of the tendencies of state variables
3016	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
3017	REAL, DIMENSION (klon), INTENT(OUT) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
3018	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
3019	REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_acol, d_adens_icol, d_adens_bnd
3020
3021
3022	!! internal variables
3023
3024	INTEGER :: i, k
3025	REAL, DIMENSION (klon) :: iwdens, isigmaw !! inactive wake density and fractional area
3026	!! REAL, DIMENSION (klon) :: d_arad, d_irad
3027	REAL, DIMENSION (klon) :: igfl !! LgI = gust front length of inactive wakes
3028	!! per unit area
3029	REAL, DIMENSION (klon) :: s_wk !! mean area of individual wakes
3030	REAL, DIMENSION (klon) :: as_wk !! mean area of individual active wakes
3031	REAL, DIMENSION (klon) :: is_wk !! mean area of individual inactive wakes
3032	REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes
3033	REAL :: tau_wk_inv_min
3034	REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes
3035	REAL :: d_wdens_targ, d_sigmaw_targ
3036
3037
3038	!! Equations
3039	!! ---------
3040	!! Gust fronts:
3041	!! Lg_A = 2 pi r_A A
3042	!! Lg_I = 2 pi r_I I
3043	!! Lg = 2 pi r D
3044	!!
3045	!! Areas:
3046	!! s = pi r^2
3047	!! s_A = pi r_A^2
3048	!! s_I = pi r_I^2
3049	!!
3050	!! Life expectancy:
3051	!! tau_I = 3 C* ((C/Ct)^3/2 - 1) / r_I
3052	!!
3053	!! Time deratives:
3054	!! dD/dt = B - (D-A)/tau_I - 2 Lg C* D
3055	!! dA/dt = B - A/tau_A + 2 Lg_I C* (D-A) - 2 Lg_A C* A
3056	!! dsigma/dt = B a0 - sigma_I/tau_I + Lg C* - 2 Lg_I C* (D-A) (2 s_I - a0)
3057	!! dsigma_A/dt = B a0 - sigma_A/tau_A + Lg_A C* + (Lg_A I + Lg_I A) C* s_I + 2 Lg_I C* I a0
3058	!!
3059
3060	DO i = 1, klon
3061	IF (wk_adv(i)) THEN
3062	iwdens(i) = wdens(i) - awdens(i)
3063	isigmaw(i) = sigmaw(i) - asigmaw(i)
3064	!
3065	arad_wk(i) = max( sqrt(asigmaw(i)/(3.14*awdens(i))) , rzero)
3066	irad_wk(i) = max( sqrt((sigmaw(i)-asigmaw(i))/ &
3067	(3.14*max(smallestreal,(wdens(i)-awdens(i))))), rzero)
3068	rad_wk(i) = (awdens(i)arad_wk(i)+(wdens(i)-awdens(i))irad_wk(i))/wdens(i)
3069	!
3070	s_wk(i) = 3.14rad_wk(i)*2
3071	as_wk(i) = 3.14arad_wk(i)*2
3072	is_wk(i) = 3.14irad_wk(i)*2
3073	!
3074	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
3075	agfl(i) = 2.sqrt(3.14awdens(i)*asigmaw(i))
3076	igfl(i) = gfl(i) - agfl(i)
3077	ENDIF
3078	ENDDO
3079
3080
3081	DO i = 1, klon
3082	IF (wk_adv(i)) THEN
3083	tau_wk_inv(i) = max( (3.cstar(i))/(irad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
3084	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
3085	tau_prime(i) = tau_cv
3086
3087	d_sig_gen(i) = wgen(i)*aa0
3088	d_sig_death(i) = - isigmaw(i)*tau_wk_inv_min
3089	d_sig_col(i) = - 2.igfl(i)cstar(i)iwdens(i)(2.*is_wk(i)-aa0)
3090	d_sig_spread(i) = gfl(i)*cstar(i)
3091	!
3092	d_sig_gen(i) = d_sig_gen(i)*dtimesub
3093	d_sig_death(i) = d_sig_death(i)*dtimesub
3094	d_sig_col(i) = d_sig_col(i)*dtimesub
3095	d_sig_spread(i) = d_sig_spread(i)*dtimesub
3096	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
3097	#ifdef IOPHYS_WK
3098	IF (phys_sub) call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)
3099	#endif
3100
3101
3102	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
3103	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
3104	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
3105	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
3106	d_sigmaw(i) = d_sigmaw_targ
3107	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
3108	#ifdef IOPHYS_WK
3109	IF (phys_sub) THEN
3110	call iophys_ecrit('tauwk_inv',1,'tau_wk_inv_min','',tau_wk_inv_min)
3111	call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw)
3112	call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen)
3113	call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death)
3114	call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col)
3115	call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread)
3116	call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd)
3117	ENDIF
3118	#endif
3119	d_asig_death(i) = - asigmaw(i)/tau_prime(i)
3120	d_asig_aicol(i) = (agfl(i)iwdens(i) + igfl(i)awdens(i))cstar(i)is_wk(i)
3121	d_asig_iicol(i) = 2.igfl(i)cstar(i)iwdens(i)aa0
3122	d_asig_spread(i) = agfl(i)*cstar(i)
3123	!
3124	d_asig_death(i) = d_asig_death(i)*dtimesub
3125	d_asig_aicol(i) = d_asig_aicol(i)*dtimesub
3126	d_asig_iicol(i) = d_asig_iicol(i)*dtimesub
3127	d_asig_spread(i) = d_asig_spread(i)*dtimesub
3128	d_asigmaw(i) = d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i)
3129	#ifdef IOPHYS_WK
3130	IF (phys_sub) call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)
3131	#endif
3132
3133	d_sigmaw_targ = min(max(d_asigmaw(i),-asigmaw(i)), sigmaw(i)-asigmaw(i))
3134	!! d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
3135	d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i)
3136	d_asigmaw(i) = d_sigmaw_targ
3137	#ifdef IOPHYS_WK
3138	IF (phys_sub) THEN
3139	call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw)
3140	call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death)
3141	call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol)
3142	call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol)
3143	call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread)
3144	call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd)
3145	ENDIF
3146	#endif
3147	d_dens_gen(i) = wgen(i)
3148	d_dens_death(i) = - iwdens(i)*tau_wk_inv_min
3149	d_dens_col(i) = - 2.gfl(i)cstar(i)*wdens(i)
3150	!
3151	d_dens_gen(i) = d_dens_gen(i)*dtimesub
3152	d_dens_death(i) = d_dens_death(i)*dtimesub
3153	d_dens_col(i) = d_dens_col(i)*dtimesub
3154	d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
3155	!!
3156	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
3157	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
3158	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
3159	d_wdens(i) = d_wdens_targ
3160	#ifdef IOPHYS_WK
3161	IF (phys_sub) THEN
3162	call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens)
3163	call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen)
3164	call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death)
3165	call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col)
3166	ENDIF
3167	#endif
3168
3169	d_adens_death(i) = -awdens(i)/tau_prime(i)
3170	d_adens_icol(i) = 2.igfl(i)cstar(i)*iwdens(i)
3171	d_adens_acol(i) = - 2.agfl(i)cstar(i)*awdens(i)
3172	!
3173	d_adens_death(i) = d_adens_death(i)*dtimesub
3174	d_adens_icol(i) = d_adens_icol(i)*dtimesub
3175	d_adens_acol(i) = d_adens_acol(i)*dtimesub
3176	d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)
3177	#ifdef IOPHYS_WK
3178	IF (phys_sub) THEN
3179	call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens)
3180	call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death)
3181	call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol)
3182	call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol)
3183	ENDIF
3184	#endif
3185	d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
3186	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
3187	d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
3188	d_awdens(i) = d_wdens_targ
3189
3190	!! d_irad(i) = (d_sigmaw(i)-d_asigmaw(i)-isigmaw(i)*(d_wdens(i)-awdens(i))/iwdens(i)) / &
3191	!! max(smallestreal,(2.3.14iwdens(i)*irad_wk(i)))
3192	!! d_arad(i) = (d_asigmaw(i)-asigmaw(i)*d_awdens(i)/awdens(i)) / &
3193	!! max(smallestreal,(2.3.14awdens(i)*arad_wk(i)))
3194	!! d_irad(i) = d_irad(i)*dtimesub
3195	!! d_arad(i) = d_arad(i)*dtimesub
3196	!! call iophys_ecrit('d_irad',1,'d_irad','m',d_irad)
3197	!! call iophys_ecrit('d_airad',1,'d_arad','m',d_arad)
3198	!!
3199	ENDIF
3200	ENDDO
3201
3202	IF (prt_level >= 10) THEN
3203	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) ', &
3204	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1)
3205	print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
3206	wdens(1), awdens(1), d_awdens(1)
3207	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
3208	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
3209	ENDIF
3210	sigmaw=sigmaw+d_sigmaw
3211	asigmaw=asigmaw+d_asigmaw
3212	wdens=wdens+d_wdens
3213	awdens=awdens+d_awdens
3214
3215	RETURN
3216	END SUBROUTINE wake_popdyn_3
3217
3218	END MODULE lmdz_wake

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

Download in other formats:

Original Format