Context Navigation

lmdz_wake.F90 @ 5051

Last change on this file since 5051 was 4908, checked in by jyg, 7 months ago
New wake height formulation in the wake parametrization
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.1 KB

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

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

Download in other formats:

Original Format