Context Navigation

lmdz_wake.F90 @ 5218

Last change on this file since 5218 was 5193, checked in by abarral, 3 months ago
Merge r5181, r5188 from trunk Lint lmdz_wake.F90, lmdz_wake_ini.F90
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.7 KB

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

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

Download in other formats:

Original Format