Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmdz_wake.F90 @ 5024

Last change on this file since 5024 was 4727, checked in by idelkadi, 13 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
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.8 KB

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

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

Download in other formats:

Original Format