Context Navigation

lmdz_wake.F90 @ 4678

Last change on this file since 4678 was 4588, checked in by fhourdin, 11 months ago
Nouveaux noms et encapsulation module pour 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: 89.2 KB

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

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

Download in other formats:

Original Format