Context Navigation

lmdz_wake.F90 @ 5212

Last change on this file since 5212 was 5203, checked in by Laurent Fairhead, 30 hours ago
Merge with trunk revision 5202 before reintegration to trunk
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 120.3 KB

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

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

Download in other formats:

Original Format