Context Navigation

wake.F90 @ 4570

Last change on this file since 4570 was 4453, checked in by fhourdin, 21 months ago
Changement d'un nom de param
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.1 KB

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

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

Download in other formats:

Original Format