Context Navigation

wake.F90 @ 4584

Last change on this file since 4584 was 4584, checked in by Laurent Fairhead, 11 months ago
Merged trunk revisions from r4443 to r4582 (HEAD) into 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.1 KB

Line
1
2	! $Id: wake.F90 4584 2023-06-23 14:10:41Z fairhead $
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