Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmdz_wake2.f90 @ 5821

Last change on this file since 5821 was 5804, checked in by fhourdin, 3 months ago
Séparation de lmdz_wake en petits fichiers (JYG&FH)
File size: 85.7 KB

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

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

Download in other formats:

Original Format