Context Navigation

lmdz_wake.F90 @ 5179

Last change on this file since 5179 was 5160, checked in by abarral, 16 months ago
Put .h into modules
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: 119.9 KB

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

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

Download in other formats:

Original Format