Context Navigation

lmdz_wake.f90 @ 5441

Last change on this file since 5441 was 5400, checked in by evignon, 12 days ago
ajout de omp threadprivate manquants
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 120.2 KB

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

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

Download in other formats:

Original Format