Context Navigation

lmdz_wake.F90 @ 4771

Last change on this file since 4771 was 4744, checked in by jyg, 12 months ago
Implementation of a two radii model for wake population dynamics. It is activated with : iflag_wk_pop_dyn=3
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 115.8 KB

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

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

Download in other formats:

Original Format