Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmdz_wake3.f90 @ 5804

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

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

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

Download in other formats:

Original Format