Context Navigation

wake.F90 @ 4294

Last change on this file since 4294 was 4294, checked in by fhourdin, 21 months ago
Nouvelle version des wake / Jean-Yves et Lamine
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: 88.9 KB

Line
1
2	! $Id: wake.F90 4294 2022-10-06 22:24:21Z fhourdin $
3
4
5	SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, &
6	tenv0, qe0, omgb, &
7	dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
8	sigd_con, Cin, &
9	deltatw, deltaqw, sigmaw, awdens, wdens, &
10	dth, hw, wape, fip, gfl, &
11	dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
12	dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
13	d_deltat_gw, & ! tendencies
14	d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2) ! tendencies
15
16
17	! **************************************************************
18	! *
19	! WAKE *
20	! retour a un Pupper fixe *
21	! *
22	! written by : GRANDPEIX Jean-Yves 09/03/2000 *
23	! modified by : ROEHRIG Romain 01/29/2007 *
24	! **************************************************************
25
26
27	USE wake_ini_mod , ONLY : wake_ini
28	USE wake_ini_mod , ONLY : prt_level,epsim1,RG,RD
29	USE wake_ini_mod , ONLY : stark, wdens_ref, coefgw, alpk, pupperbyphs
30	USE wake_ini_mod , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
31	USE wake_ini_mod , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin
32	USE wake_ini_mod , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
33	USE wake_ini_mod , ONLY : iflag_wk_profile
34
35
36	IMPLICIT NONE
37	! ============================================================================
38
39
40	! But : Decrire le comportement des poches froides apparaissant dans les
41	! grands systemes convectifs, et fournir l'energie disponible pour
42	! le declenchement de nouvelles colonnes convectives.
43
44	! State variables :
45	! deltatw : temperature difference between wake and off-wake regions
46	! deltaqw : specific humidity difference between wake and off-wake regions
47	! sigmaw : fractional area covered by wakes.
48	! wdens : number of wakes per unit area
49
50	! Variable de sortie :
51
52	! wape : WAke Potential Energy
53	! fip : Front Incident Power (W/m2) - ALP
54	! gfl : Gust Front Length per unit area (m-1)
55	! dtls : large scale temperature tendency due to wake
56	! dqls : large scale humidity tendency due to wake
57	! hw : wake top hight (given by hw*deltatw(1)/2=wape)
58	! dp_omgb : vertical gradient of large scale omega
59	! awdens : densite de poches actives
60	! wdens : densite de poches
61	! omgbdth: flux of Delta_Theta transported by LS omega
62	! dtKE : differential heating (wake - unpertubed)
63	! dqKE : differential moistening (wake - unpertubed)
64	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
65	! dp_deltomg : vertical gradient of omg (s-1)
66	! wkspread : spreading term in d_t_wake and d_q_wake
67	! deltatw : updated temperature difference (T_w-T_u).
68	! deltaqw : updated humidity difference (q_w-q_u).
69	! sigmaw : updated wake fractional area.
70	! d_deltat_gw : delta T tendency due to GW
71
72	! Variables d'entree :
73
74	! aire : aire de la maille
75	! tenv0 : temperature dans l'environnement (K)
76	! qe0 : humidite dans l'environnement (kg/kg)
77	! omgb : vitesse verticale moyenne sur la maille (Pa/s)
78	! dtdwn: source de chaleur due aux descentes (K/s)
79	! dqdwn: source d'humidite due aux descentes (kg/kg/s)
80	! dta : source de chaleur due courants satures et detrain (K/s)
81	! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
82	! wgen : number of wakes generated per unit area and per sec (/m^2/s)
83	! amdwn: flux de masse total des descentes, par unite de
84	! surface de la maille (kg/m2/s)
85	! amup : flux de masse total des ascendances, par unite de
86	! surface de la maille (kg/m2/s)
87	! sigd_con:
88	! Cin : convective inhibition
89	! p : pressions aux milieux des couches (Pa)
90	! ph : pressions aux interfaces (Pa)
91	! pi : (p/p_0)**kapa (adim)
92	! dtime: increment temporel (s)
93
94	! Variables internes :
95
96	! rhow : masse volumique de la poche froide
97	! rho : environment density at P levels
98	! rhoh : environment density at Ph levels
99	! tenv : environment temperature \| may change within
100	! qe : environment humidity \| sub-time-stepping
101	! the : environment potential temperature
102	! thu : potential temperature in undisturbed area
103	! tu : temperature in undisturbed area
104	! qu : humidity in undisturbed area
105	! dp_omgb: vertical gradient og LS omega
106	! omgbw : wake average vertical omega
107	! dp_omgbw: vertical gradient of omgbw
108	! omgbdq : flux of Delta_q transported by LS omega
109	! dth : potential temperature diff. wake-undist.
110	! th1 : first pot. temp. for vertical advection (=thu)
111	! th2 : second pot. temp. for vertical advection (=thw)
112	! q1 : first humidity for vertical advection
113	! q2 : second humidity for vertical advection
114	! d_deltatw : terme de redistribution pour deltatw
115	! d_deltaqw : terme de redistribution pour deltaqw
116	! deltatw0 : deltatw initial
117	! deltaqw0 : deltaqw initial
118	! hw0 : wake top hight (defined as the altitude at which deltatw=0)
119	! amflux : horizontal mass flux through wake boundary
120	! wdens_ref: initial number of wakes per unit area (3D) or per
121	! unit length (2D), at the beginning of each time step
122	! Tgw : 1 sur la periode de onde de gravite
123	! Cgw : vitesse de propagation de onde de gravite
124	! LL : distance entre 2 poches
125
126	! -------------------------------------------------------------------------
127	! Declaration de variables
128	! -------------------------------------------------------------------------
129
130
131	! Arguments en entree
132	! --------------------
133
134	INTEGER, INTENT(IN) :: klon,klev
135	INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf
136	REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi
137	REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph
138	REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb
139	REAL, INTENT(IN) :: dtime
140	REAL, DIMENSION (klon, klev), INTENT(IN) :: tenv0, qe0
141	REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn
142	REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup
143	REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa
144	REAL, DIMENSION (klon), INTENT(IN) :: wgen
145	REAL, DIMENSION (klon), INTENT(IN) :: sigd_con
146	REAL, DIMENSION (klon), INTENT(IN) :: Cin
147
148	!
149	! Input/Output
150	! State variables
151	REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw
152	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw
153	REAL, DIMENSION (klon), INTENT(INOUT) :: awdens
154	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens
155
156	! Sorties
157	! --------
158
159	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth
160	REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu
161	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls
162	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke
163	REAL, DIMENSION (klon, klev), INTENT(OUT) :: wkspread ! unused (jyg)
164	REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg
165	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg
166	REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar
167	INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw
168	! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
169	! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
170	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw
171	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2
172	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_awdens2, d_wdens2
173
174	! Variables internes
175	! -------------------
176
177	! Variables a fixer
178
179	REAL :: delta_t_min
180	INTEGER :: nsub
181	REAL :: dtimesub
182	REAL :: wdens0
183	! IM 080208
184	LOGICAL, DIMENSION (klon) :: gwake
185
186	! Variables de sauvegarde
187	REAL, DIMENSION (klon, klev) :: deltatw0
188	REAL, DIMENSION (klon, klev) :: deltaqw0
189	REAL, DIMENSION (klon, klev) :: tenv, qe
190	!! REAL, DIMENSION (klon) :: sigmaw1
191
192	! Variables liees a la dynamique de population 1
193	REAL, DIMENSION(klon) :: act
194	REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv
195	REAL, DIMENSION(klon) :: f_shear
196	REAL, DIMENSION(klon) :: drdt
197
198	! Variables liees a la dynamique de population 2
199	REAL, DIMENSION(klon) :: cont_fact
200
201	!! REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col
202	REAL, DIMENSION(klon) :: wape1_act, wape2_act
203	LOGICAL, DIMENSION (klon) :: kill_wake
204	REAL :: drdt_pos
205	REAL :: tau_wk_inv_min
206	! Some components of the tendencies of state variables
207	REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
208	REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
209	REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
210
211	! Variables pour les GW
212	REAL, DIMENSION (klon) :: ll
213	REAL, DIMENSION (klon, klev) :: n2
214	REAL, DIMENSION (klon, klev) :: cgw
215	REAL, DIMENSION (klon, klev) :: tgw
216
217	! Variables liees au calcul de hw
218	REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
219	REAL, DIMENSION (klon) :: sum_dth
220	REAL, DIMENSION (klon) :: dthmin
221	REAL, DIMENSION (klon) :: z, dz, hw0
222	INTEGER, DIMENSION (klon) :: ktop, kupper
223
224	! Variables liees au test de la forme triangulaire du profil de Delta_theta
225	REAL, DIMENSION (klon) :: sum_half_dth
226	REAL, DIMENSION (klon) :: dz_half
227
228	! Sub-timestep tendencies and related variables
229	REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw
230	REAL, DIMENSION (klon, klev) :: d_tenv, d_qe
231	REAL, DIMENSION (klon) :: d_awdens, d_wdens, d_sigmaw
232	REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
233	REAL, DIMENSION (klon) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
234	REAL, DIMENSION (klon) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
235	REAL, DIMENSION (klon) :: alpha, alpha_tot
236	REAL, DIMENSION (klon) :: q0_min, q1_min
237	LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw
238
239	! Autres variables internes
240	INTEGER ::isubstep, k, i, igout
241
242	REAL :: sigmaw_targ
243	REAL :: wdens_targ
244	REAL :: d_sigmaw_targ
245	REAL :: d_wdens_targ
246
247	REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
248	REAL, DIMENSION (klon) :: sum_dq, sum_rho
249	REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
250	REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
251	REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
252	REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
253
254	REAL, DIMENSION (klon, klev) :: rho, rhow
255	REAL, DIMENSION (klon, klev+1) :: rhoh
256	REAL, DIMENSION (klon, klev) :: rhow_moyen
257	REAL, DIMENSION (klon, klev) :: zh
258	REAL, DIMENSION (klon, klev+1) :: zhh
259	REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
260
261	REAL, DIMENSION (klon, klev) :: the, thu
262
263	REAL, DIMENSION (klon, klev) :: omgbw
264	REAL, DIMENSION (klon) :: pupper
265	REAL, DIMENSION (klon) :: omgtop
266	REAL, DIMENSION (klon, klev) :: dp_omgbw
267	REAL, DIMENSION (klon) :: ztop, dztop
268	REAL, DIMENSION (klon, klev) :: alpha_up
269
270	REAL, DIMENSION (klon) :: rre1, rre2
271	REAL :: rrd1, rrd2
272	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
273	REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
274	REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
275	REAL, DIMENSION (klon, klev) :: omgbdq
276
277	REAL, DIMENSION (klon) :: ff, gg
278	REAL, DIMENSION (klon) :: wape2, cstar2, heff
279
280	REAL, DIMENSION (klon, klev) :: crep
281
282	REAL, DIMENSION (klon, klev) :: ppi
283
284	! cc nrlmd
285	REAL, DIMENSION (klon) :: death_rate
286	!! REAL, DIMENSION (klon) :: nat_rate
287	REAL, DIMENSION (klon, klev) :: entr
288	REAL, DIMENSION (klon, klev) :: detr
289
290	REAL, DIMENSION(klon) :: sigmaw_in ! pour les prints
291	REAL, DIMENSION(klon) :: awdens_in, wdens_in ! pour les prints
292
293	write(81) klon,klev,znatsurf,p,pi,ph,omgb,dtime,tenv0,qe0,dtdwn,dqdwn,amdwn,amup,dta,dqa,wgen,sigd_con,cin,deltatw,deltaqw,sigmaw,awdens,wdens
294	! -------------------------------------------------------------------------
295	! Initialisations
296	! -------------------------------------------------------------------------
297	! ALON = 3.e5
298	! alon = 1.E6
299
300	! Provisionnal; to be suppressed when f_shear is parameterized
301	f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear
302
303
304	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
305
306	! coefgw : Coefficient pour les ondes de gravite
307	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
308	! wdens : Densite surfacique de poche froide
309	! -------------------------------------------------------------------------
310
311	! cc nrlmd coefgw=10
312	! coefgw=1
313	! wdens0 = 1.0/(alon**2)
314	! cc nrlmd wdens = 1.0/(alon**2)
315	! cc nrlmd stark = 0.50
316	! CRtest
317	! cc nrlmd alpk=0.1
318	! alpk = 1.0
319	! alpk = 0.5
320	! alpk = 0.05
321	print *,'XXXX dtime input ', dtime
322	igout = klon/2+1/klon
323
324	IF (iflag_wk_pop_dyn == 0) THEN
325	! Initialisation de toutes des densites a wdens_ref.
326	! Les densites peuvent evoluer si les poches debordent
327	! (voir au tout debut de la boucle sur les substeps)
328	!jyg<
329	!! wdens(:) = wdens_ref
330	DO i = 1,klon
331	wdens(i) = wdens_ref(znatsurf(i)+1)
332	ENDDO
333	!>jyg
334	ENDIF ! (iflag_wk_pop_dyn == 0)
335
336	! print*,'stark',stark
337	! print*,'alpk',alpk
338	! print*,'wdens',wdens
339	! print*,'coefgw',coefgw
340	! cc
341	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
342	! -------------------------------------------------------------------------
343
344	delta_t_min = 0.2
345
346	! 1. - Save initial values, initialize tendencies, initialize output fields
347	! ------------------------------------------------------------------------
348
349	!jyg<
350	!! DO k = 1, klev
351	!! DO i = 1, klon
352	!! ppi(i, k) = pi(i, k)
353	!! deltatw0(i, k) = deltatw(i, k)
354	!! deltaqw0(i, k) = deltaqw(i, k)
355	!! tenv(i, k) = tenv0(i, k)
356	!! qe(i, k) = qe0(i, k)
357	!! dtls(i, k) = 0.
358	!! dqls(i, k) = 0.
359	!! d_deltat_gw(i, k) = 0.
360	!! d_tenv(i, k) = 0.
361	!! d_qe(i, k) = 0.
362	!! d_deltatw(i, k) = 0.
363	!! d_deltaqw(i, k) = 0.
364	!! ! IM 060508 beg
365	!! d_deltatw2(i, k) = 0.
366	!! d_deltaqw2(i, k) = 0.
367	!! ! IM 060508 end
368	!! END DO
369	!! END DO
370	ppi(:,:) = pi(:,:)
371	deltatw0(:,:) = deltatw(:,:)
372	deltaqw0(:,:) = deltaqw(:,:)
373	tenv(:,:) = tenv0(:,:)
374	qe(:,:) = qe0(:,:)
375	dtls(:,:) = 0.
376	dqls(:,:) = 0.
377	d_deltat_gw(:,:) = 0.
378	d_tenv(:,:) = 0.
379	d_qe(:,:) = 0.
380	d_deltatw(:,:) = 0.
381	d_deltaqw(:,:) = 0.
382	d_deltatw2(:,:) = 0.
383	d_deltaqw2(:,:) = 0.
384
385	IF (iflag_wk_act == 0) THEN
386	act(:) = 0.
387	ELSEIF (iflag_wk_act == 1) THEN
388	act(:) = 1.
389	ENDIF
390
391	!! DO i = 1, klon
392	!! sigmaw_in(i) = sigmaw(i)
393	!! END DO
394	sigmaw_in(:) = sigmaw(:)
395	!>jyg
396	!
397	IF (iflag_wk_pop_dyn >= 1) THEN
398	awdens_in(:) = awdens(:)
399	wdens_in(:) = wdens(:)
400	!! wdens(:) = wdens(:) + wgen(:)*dtime
401	!! d_wdens2(:) = wgen(:)*dtime
402	!! ELSE
403	ENDIF ! (iflag_wk_pop_dyn >= 1)
404
405
406	! sigmaw1=sigmaw
407	! IF (sigd_con.GT.sigmaw1) THEN
408	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
409	! ENDIF
410	IF (iflag_wk_pop_dyn >= 1) THEN
411	DO i = 1, klon
412	d_dens_gen2(i) = 0.
413	d_dens_death2(i) = 0.
414	d_dens_col2(i) = 0.
415	d_awdens2(i) = 0.
416	!
417	wdens_targ = max(wdens(i),wdensmin)
418	d_dens_bnd2(i) = wdens_targ - wdens(i)
419	d_wdens2(i) = wdens_targ - wdens(i)
420	wdens(i) = wdens_targ
421	END DO
422	IF (iflag_wk_pop_dyn == 2) THEN
423	DO i = 1, klon
424	d_adens_death2(i) = 0.
425	d_adens_icol2(i) = 0.
426	d_adens_acol2(i) = 0.
427	!
428	wdens_targ = min(max(awdens(i),0.),wdens(i))
429	d_adens_bnd2(i) = wdens_targ - awdens(i)
430	d_awdens2(i) = wdens_targ - awdens(i)
431	awdens(i) = wdens_targ
432	END DO
433	ENDIF ! (iflag_wk_pop_dyn == 2)
434	ELSE
435	DO i = 1, klon
436	d_awdens2(i) = 0.
437	d_wdens2(i) = 0.
438	END DO
439	ENDIF ! (iflag_wk_pop_dyn >= 1)
440	!
441	DO i = 1, klon
442	! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
443	!jyg<
444	!! sigmaw(i) = amax1(sigmaw(i), sigmad)
445	!! sigmaw(i) = amin1(sigmaw(i), 0.99)
446	d_sig_gen2(i) = 0.
447	d_sig_death2(i) = 0.
448	d_sig_col2(i) = 0.
449	d_sig_spread2(i)= 0.
450	sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
451	d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
452	d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
453	print *,'XXXX1 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
454	sigmaw(i) = sigmaw_targ
455	!>jyg
456	END DO
457
458	wape(:) = 0.
459	wape2(:) = 0.
460	d_sigmaw(:) = 0.
461	ktopw(:) = 0
462	!
463	!<jyg
464	dth(:,:) = 0.
465	tu(:,:) = 0.
466	qu(:,:) = 0.
467	dtke(:,:) = 0.
468	dqke(:,:) = 0.
469	wkspread(:,:) = 0.
470	omgbdth(:,:) = 0.
471	omg(:,:) = 0.
472	dp_omgb(:,:) = 0.
473	dp_deltomg(:,:) = 0.
474	hw(:) = 0.
475	wape(:) = 0.
476	fip(:) = 0.
477	gfl(:) = 0.
478	cstar(:) = 0.
479	ktopw(:) = 0
480	!
481	! Vertical advection local variables
482	omgbw(:,:) = 0.
483	omgtop(:) = 0
484	dp_omgbw(:,:) = 0.
485	omgbdq(:,:) = 0.
486
487	!>jyg
488	!
489	IF (prt_level>=10) THEN
490	PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
491	PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
492	PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
493	PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
494	PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
495	PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
496	PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
497	PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
498	PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
499	ENDIF
500
501	! 2. - Prognostic part
502	! --------------------
503
504
505	! 2.1 - Undisturbed area and Wake integrals
506	! ---------------------------------------------------------
507
508	DO i = 1, klon
509	z(i) = 0.
510	ktop(i) = 0
511	kupper(i) = 0
512	sum_thu(i) = 0.
513	sum_tu(i) = 0.
514	sum_qu(i) = 0.
515	sum_thvu(i) = 0.
516	sum_dth(i) = 0.
517	sum_dq(i) = 0.
518	sum_rho(i) = 0.
519	sum_dtdwn(i) = 0.
520	sum_dqdwn(i) = 0.
521
522	av_thu(i) = 0.
523	av_tu(i) = 0.
524	av_qu(i) = 0.
525	av_thvu(i) = 0.
526	av_dth(i) = 0.
527	av_dq(i) = 0.
528	av_rho(i) = 0.
529	av_dtdwn(i) = 0.
530	av_dqdwn(i) = 0.
531	END DO
532
533	! Distance between wakes
534	DO i = 1, klon
535	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
536	END DO
537	! Potential temperatures and humidity
538	! ----------------------------------------------------------
539	DO k = 1, klev
540	DO i = 1, klon
541	! write(,)'wake 1',i,k,RD,tenv(i,k)
542	rho(i, k) = p(i, k)/(RD*tenv(i,k))
543	! write(,)'wake 2',rho(i,k)
544	IF (k==1) THEN
545	! write(,)'wake 3',i,k,rd,tenv(i,k)
546	rhoh(i, k) = ph(i, k)/(RD*tenv(i,k))
547	! write(,)'wake 4',i,k,rd,tenv(i,k)
548	zhh(i, k) = 0
549	ELSE
550	! write(,)'wake 5',rd,(tenv(i,k)+tenv(i,k-1))
551	rhoh(i, k) = ph(i, k)2./(RD(tenv(i,k)+tenv(i,k-1)))
552	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
553	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
554	END IF
555	! write(,)'wake 7',ppi(i,k)
556	the(i, k) = tenv(i, k)/ppi(i, k)
557	thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
558	tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i)
559	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
560	! write(,)'wake 8',(RD*(tenv(i,k)+deltatw(i,k)))
561	rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k)))
562	dth(i, k) = deltatw(i, k)/ppi(i, k)
563	END DO
564	END DO
565
566	DO k = 1, klev - 1
567	DO i = 1, klon
568	IF (k==1) THEN
569	n2(i, k) = 0
570	ELSE
571	n2(i, k) = amax1(0., -RG*2/the(i,k)rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
572	(p(i,k+1)-p(i,k-1)))
573	END IF
574	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
575
576	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
577	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
578	END DO
579	END DO
580
581	DO i = 1, klon
582	n2(i, klev) = 0
583	zh(i, klev) = 0
584	cgw(i, klev) = 0
585	tgw(i, klev) = 0
586	END DO
587
588	! Calcul de la masse volumique moyenne de la colonne (bdlmd)
589	! -----------------------------------------------------------------
590
591	DO k = 1, klev
592	DO i = 1, klon
593	epaisseur1(i, k) = 0.
594	epaisseur2(i, k) = 0.
595	END DO
596	END DO
597
598	DO i = 1, klon
599	epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.
600	epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.
601	rhow_moyen(i, 1) = rhow(i, 1)
602	END DO
603
604	DO k = 2, klev
605	DO i = 1, klon
606	epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG) + 1.
607	epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
608	rhow_moyen(i, k) = (rhow_moyen(i,k-1)epaisseur2(i,k-1)+rhow(i,k) &
609	epaisseur1(i,k))/epaisseur2(i, k)
610	END DO
611	END DO
612
613
614	! Choose an integration bound well above wake top
615	! -----------------------------------------------------------------
616
617	! Determine Wake top pressure (Ptop) from buoyancy integral
618	! --------------------------------------------------------
619
620	! -1/ Pressure of the level where dth becomes less than delta_t_min.
621
622	DO i = 1, klon
623	ptop_provis(i) = ph(i, 1)
624	END DO
625	DO k = 2, klev
626	DO i = 1, klon
627
628	! IM v3JYG; ptop_provis(i).LT. ph(i,1)
629
630	IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
631	ptop_provis(i)==ph(i,1)) THEN
632	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &
633	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
634	END IF
635	END DO
636	END DO
637
638	! -2/ dth integral
639
640	DO i = 1, klon
641	sum_dth(i) = 0.
642	dthmin(i) = -delta_t_min
643	z(i) = 0.
644	END DO
645
646	DO k = 1, klev
647	DO i = 1, klon
648	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
649	IF (dz(i)>0) THEN
650	z(i) = z(i) + dz(i)
651	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
652	dthmin(i) = amin1(dthmin(i), dth(i,k))
653	END IF
654	END DO
655	END DO
656
657	! -3/ height of triangle with area= sum_dth and base = dthmin
658
659	DO i = 1, klon
660	hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
661	hw0(i) = amax1(hwmin, hw0(i))
662	END DO
663
664	! -4/ now, get Ptop
665
666	DO i = 1, klon
667	z(i) = 0.
668	ptop(i) = ph(i, 1)
669	END DO
670
671	DO k = 1, klev
672	DO i = 1, klon
673	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw0(i)-z(i))
674	IF (dz(i)>0) THEN
675	z(i) = z(i) + dz(i)
676	ptop(i) = ph(i, k) - rho(i, k)RGdz(i)
677	END IF
678	END DO
679	END DO
680
681	IF (prt_level>=10) THEN
682	PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
683	ENDIF
684
685
686	! -5/ Determination de ktop et kupper
687
688	CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
689
690	DO k = klev, 1, -1
691	DO i = 1, klon
692	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
693	END DO
694	END DO
695	!print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
696
697
698
699	! -6/ Correct ktop and ptop
700
701	DO i = 1, klon
702	ptop_new(i) = ptop(i)
703	END DO
704	DO k = klev, 2, -1
705	DO i = 1, klon
706	IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
707	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
708	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
709	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
710	END IF
711	END DO
712	END DO
713
714	DO i = 1, klon
715	ptop(i) = ptop_new(i)
716	END DO
717
718	DO k = klev, 1, -1
719	DO i = 1, klon
720	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
721	END DO
722	END DO
723
724	IF (prt_level>=10) THEN
725	PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
726	ENDIF
727
728	! -5/ Set deltatw & deltaqw to 0 above kupper
729
730	DO k = 1, klev
731	DO i = 1, klon
732	IF (k>=kupper(i)) THEN
733	deltatw(i, k) = 0.
734	deltaqw(i, k) = 0.
735	d_deltatw2(i,k) = -deltatw0(i,k)
736	d_deltaqw2(i,k) = -deltaqw0(i,k)
737	END IF
738	END DO
739	END DO
740
741
742	! Vertical gradient of LS omega
743
744	DO k = 1, klev
745	DO i = 1, klon
746	IF (k<=kupper(i)) THEN
747	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
748	END IF
749	END DO
750	END DO
751
752	! Integrals (and wake top level number)
753	! --------------------------------------
754
755	! Initialize sum_thvu to 1st level virt. pot. temp.
756
757	DO i = 1, klon
758	z(i) = 1.
759	dz(i) = 1.
760	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
761	sum_dth(i) = 0.
762	END DO
763
764	DO k = 1, klev
765	DO i = 1, klon
766	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
767	IF (dz(i)>0) THEN
768	z(i) = z(i) + dz(i)
769	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
770	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
771	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
772	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
773	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
774	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
775	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
776	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
777	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
778	END IF
779	END DO
780	END DO
781
782	DO i = 1, klon
783	hw0(i) = z(i)
784	END DO
785
786
787	! 2.1 - WAPE and mean forcing computation
788	! ---------------------------------------
789
790	! ---------------------------------------
791
792	! Means
793
794	DO i = 1, klon
795	av_thu(i) = sum_thu(i)/hw0(i)
796	av_tu(i) = sum_tu(i)/hw0(i)
797	av_qu(i) = sum_qu(i)/hw0(i)
798	av_thvu(i) = sum_thvu(i)/hw0(i)
799	! av_thve = sum_thve/hw0
800	av_dth(i) = sum_dth(i)/hw0(i)
801	av_dq(i) = sum_dq(i)/hw0(i)
802	av_rho(i) = sum_rho(i)/hw0(i)
803	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
804	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
805
806	wape(i) = -RGhw0(i)(av_dth(i)+ &
807	epsim1(av_thu(i)av_dq(i)+av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
808
809	END DO
810
811	! 2.2 Prognostic variable update
812	! ------------------------------
813
814	! Filter out bad wakes
815
816	DO k = 1, klev
817	DO i = 1, klon
818	IF (wape(i)<0.) THEN
819	deltatw(i, k) = 0.
820	deltaqw(i, k) = 0.
821	dth(i, k) = 0.
822	d_deltatw2(i,k) = -deltatw0(i,k)
823	d_deltaqw2(i,k) = -deltaqw0(i,k)
824	END IF
825	END DO
826	END DO
827
828	DO i = 1, klon
829	IF (wape(i)<0.) THEN
830	wape(i) = 0.
831	cstar(i) = 0.
832	hw(i) = hwmin
833	!jyg<
834	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
835	sigmaw_targ = max(sigmad, sigd_con(i))
836	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
837	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
838	print *,'XXXX2 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
839	sigmaw(i) = sigmaw_targ
840	!>jyg
841	fip(i) = 0.
842	gwake(i) = .FALSE.
843	ELSE
844	hw(i) = hw0(i)
845	cstar(i) = starksqrt(2.wape(i))
846	gwake(i) = .TRUE.
847	END IF
848	END DO
849
850
851	! Check qx and qw positivity
852	! --------------------------
853	DO i = 1, klon
854	q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), &
855	(qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
856	END DO
857	DO k = 2, klev
858	DO i = 1, klon
859	q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
860	(qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
861	IF (q1_min(i)<=q0_min(i)) THEN
862	q0_min(i) = q1_min(i)
863	END IF
864	END DO
865	END DO
866
867	DO i = 1, klon
868	ok_qx_qw(i) = q0_min(i) >= 0.
869	alpha(i) = 1.
870	alpha_tot(i) = 1.
871	END DO
872
873	IF (prt_level>=10) THEN
874	PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
875	sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
876	ENDIF
877
878
879	! C -----------------------------------------------------------------
880	! Sub-time-stepping
881	! -----------------
882
883	nsub = 10
884	dtimesub = dtime/nsub
885
886
887
888	! ------------------------------------------------------------
889	DO isubstep = 1, nsub
890	! ------------------------------------------------------------
891	CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
892
893	!print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
894
895	! wk_adv is the logical flag enabling wake evolution in the time advance
896	! loop
897	DO i = 1, klon
898	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
899	END DO
900	IF (prt_level>=10) THEN
901	PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
902	isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
903
904	ENDIF
905
906	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
907	! negatif de ktop a kupper --------
908	! cc On calcule pour cela une densite wdens0 pour laquelle on
909	! aurait un entrainement nul ---
910	!jyg<
911	! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
912	! les poches sont insuffisantes pour accueillir tout le flux de masse
913	! des descentes unsaturees. Nous faisons alors l'hypothese que la
914	! convection profonde cree directement de nouvelles poches, sans passer
915	! par les thermiques. La nouvelle valeur de wdens est alors imposee.
916
917	DO i = 1, klon
918	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
919	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
920	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
921	IF ( iflag_wk_profile == 0 ) THEN
922	omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
923	RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
924	ELSE
925	omg(i, kupper(i)+1)=0.
926	ENDIF
927	wdens0 = (sigmaw(i)/(4.3.14)) &
928	((1.-sigmaw(i))omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))cstar(i)))**(2)
929	IF (prt_level >= 10) THEN
930	print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
931	omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
932	ENDIF
933	IF (wdens(i)<=wdens0*1.1) THEN
934	IF (iflag_wk_pop_dyn >= 1) THEN
935	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
936	d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
937	ENDIF
938	wdens(i) = wdens0
939	END IF
940	END IF
941	END DO
942
943	DO i = 1, klon
944	IF (wk_adv(i)) THEN
945	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
946	rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
947	!jyg<
948	!! sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
949	sigmaw_targ = min(sigmaw(i), sigmaw_max)
950	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
951	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
952	print *,'XXXX3 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
953	sigmaw(i) = sigmaw_targ
954	!>jyg
955	END IF
956	END DO
957
958	IF (iflag_wk_pop_dyn == 1) THEN
959
960	CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
961	dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
962	d_awdens, d_wdens, d_sigmaw, &
963	iflag_wk_act, wk_adv, cin, wape, &
964	drdt, &
965	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
966	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
967	d_wdens_targ, d_sigmaw_targ)
968
969	! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
970	! Here, it has to be set to zero.
971	death_rate(:) = 0.
972
973	ELSEIF (iflag_wk_pop_dyn == 2) THEN
974	CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
975	sigmaw, wdens, awdens, & !! states variables
976	gfl, cstar, cin, wape, rad_wk, &
977	d_sigmaw, d_wdens, d_awdens, & !! tendences
978	cont_fact, &
979	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
980	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
981	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
982	death_rate(:) = 0.
983
984	ELSEIF (iflag_wk_pop_dyn == 0.) THEN
985
986	! cc nrlmd
987
988	DO i = 1, klon
989	IF (wk_adv(i)) THEN
990	! cc nrlmd Introduction du taux de mortalite des poches et
991	! test sur sigmaw_max=0.4
992	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
993	IF (sigmaw(i)>=sigmaw_max) THEN
994	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
995	ELSE
996	death_rate(i) = 0.
997	END IF
998
999	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
1000	dtimesub
1001	! $ - nat_rate(i)sigmaw(i)dtimesub
1002	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1003	! c $ death_rate(i),ktop(i),kupper(i)',
1004	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1005	! c $ death_rate(i),ktop(i),kupper(i)
1006
1007	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
1008	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
1009	! wdens = wdens0/(10.*sigmaw)
1010	! sigmaw =max(sigmaw,sigd_con)
1011	! sigmaw =max(sigmaw,sigmad)
1012	END IF
1013	END DO
1014
1015	ENDIF ! (iflag_wk_pop_dyn == 1)
1016
1017
1018	! calcul de la difference de vitesse verticale poche - zone non perturbee
1019	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
1020	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
1021	! IM 060208 au niveau k=1...
1022	!JYG 161013 Correction : maintenant omg est dimensionne a klev.
1023	DO k = 1, klev
1024	DO i = 1, klon
1025	IF (wk_adv(i)) THEN !!! nrlmd
1026	dp_deltomg(i, k) = 0.
1027	END IF
1028	END DO
1029	END DO
1030	DO k = 1, klev
1031	DO i = 1, klon
1032	IF (wk_adv(i)) THEN !!! nrlmd
1033	omg(i, k) = 0.
1034	END IF
1035	END DO
1036	END DO
1037
1038	DO i = 1, klon
1039	IF (wk_adv(i)) THEN
1040	z(i) = 0.
1041	omg(i, 1) = 0.
1042	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
1043	END IF
1044	END DO
1045
1046	DO k = 2, klev
1047	DO i = 1, klon
1048	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1049	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
1050	z(i) = z(i) + dz(i)
1051	dp_deltomg(i, k) = dp_deltomg(i, 1)
1052	omg(i, k) = dp_deltomg(i, 1)*z(i)
1053	END IF
1054	END DO
1055	END DO
1056
1057	DO i = 1, klon
1058	IF (wk_adv(i)) THEN
1059	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
1060	ztop(i) = z(i) + dztop(i)
1061	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
1062	END IF
1063	END DO
1064
1065	IF (prt_level>=10) THEN
1066	PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
1067	PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
1068	omgtop(igout), ptop(igout), ktop(igout)
1069	ENDIF
1070
1071	! -----------------
1072	! From m/s to Pa/s
1073	! -----------------
1074
1075	DO i = 1, klon
1076	IF (wk_adv(i)) THEN
1077	omgtop(i) = -rho(i, ktop(i))RGomgtop(i)
1078	dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
1079	END IF
1080	END DO
1081
1082	DO k = 1, klev
1083	DO i = 1, klon
1084	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1085	omg(i, k) = -rho(i, k)RGomg(i, k)
1086	dp_deltomg(i, k) = dp_deltomg(i, 1)
1087	END IF
1088	END DO
1089	END DO
1090
1091	! raccordement lineaire de omg de ptop a pupper
1092
1093	DO i = 1, klon
1094	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
1095	IF ( iflag_wk_profile == 0 ) THEN
1096	omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
1097	RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
1098	ELSE
1099	omg(i, kupper(i)+1) = 0.
1100	ENDIF
1101	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
1102	(ptop(i)-pupper(i))
1103	END IF
1104	END DO
1105
1106	! c DO i=1,klon
1107	! c print*,'Pente entre 0 et kupper (reference)'
1108	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
1109	! c print*,'Pente entre ktop et kupper'
1110	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
1111	! c ENDDO
1112	! c
1113	DO k = 1, klev
1114	DO i = 1, klon
1115	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
1116	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
1117	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
1118	END IF
1119	END DO
1120	END DO
1121	!! print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
1122	! cc nrlmd
1123	! c DO i=1,klon
1124	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
1125	! c END DO
1126	! cc
1127
1128
1129	! -- Compute wake average vertical velocity omgbw
1130
1131
1132	DO k = 1, klev
1133	DO i = 1, klon
1134	IF (wk_adv(i)) THEN
1135	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
1136	END IF
1137	END DO
1138	END DO
1139	! -- and its vertical gradient dp_omgbw
1140
1141	DO k = 1, klev-1
1142	DO i = 1, klon
1143	IF (wk_adv(i)) THEN
1144	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
1145	END IF
1146	END DO
1147	END DO
1148	DO i = 1, klon
1149	IF (wk_adv(i)) THEN
1150	dp_omgbw(i, klev) = 0.
1151	END IF
1152	END DO
1153
1154	! -- Upstream coefficients for omgb velocity
1155	! -- (alpha_up(k) is the coefficient of the value at level k)
1156	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
1157	DO k = 1, klev
1158	DO i = 1, klon
1159	IF (wk_adv(i)) THEN
1160	alpha_up(i, k) = 0.
1161	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
1162	END IF
1163	END DO
1164	END DO
1165
1166	! Matrix expressing [The,deltatw] from [Th1,Th2]
1167
1168	DO i = 1, klon
1169	IF (wk_adv(i)) THEN
1170	rre1(i) = 1. - sigmaw(i)
1171	rre2(i) = sigmaw(i)
1172	END IF
1173	END DO
1174	rrd1 = -1.
1175	rrd2 = 1.
1176
1177	! -- Get [Th1,Th2], dth and [q1,q2]
1178
1179	DO k = 1, klev
1180	DO i = 1, klon
1181	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1182	dth(i, k) = deltatw(i, k)/ppi(i, k)
1183	! print *, 'VVVVwake k, the(i,k), dth(i,k), sigmaw(i) ', k, the(i,k), dth(i,k), sigmaw(i)
1184	th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
1185	th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
1186	q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
1187	q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
1188	END IF
1189	END DO
1190	END DO
1191
1192	DO i = 1, klon
1193	IF (wk_adv(i)) THEN !!! nrlmd
1194	d_th1(i, 1) = 0.
1195	d_th2(i, 1) = 0.
1196	d_dth(i, 1) = 0.
1197	d_q1(i, 1) = 0.
1198	d_q2(i, 1) = 0.
1199	d_dq(i, 1) = 0.
1200	END IF
1201	END DO
1202
1203	DO k = 2, klev
1204	DO i = 1, klon
1205	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1206	d_th1(i, k) = th1(i, k-1) - th1(i, k)
1207	d_th2(i, k) = th2(i, k-1) - th2(i, k)
1208	d_dth(i, k) = dth(i, k-1) - dth(i, k)
1209	d_q1(i, k) = q1(i, k-1) - q1(i, k)
1210	d_q2(i, k) = q2(i, k-1) - q2(i, k)
1211	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
1212	END IF
1213	END DO
1214	END DO
1215
1216	DO i = 1, klon
1217	IF (wk_adv(i)) THEN
1218	omgbdth(i, 1) = 0.
1219	omgbdq(i, 1) = 0.
1220	END IF
1221	END DO
1222
1223	DO k = 2, klev
1224	DO i = 1, klon
1225	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
1226	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
1227	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
1228	END IF
1229	END DO
1230	END DO
1231
1232	!! IF (prt_level>=10) THEN
1233	IF (prt_level>=10 .and. wk_adv(igout)) THEN
1234	PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
1235	PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
1236	PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
1237	PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
1238	ENDIF
1239
1240	! -----------------------------------------------------------------
1241	DO k = 1, klev-1
1242	DO i = 1, klon
1243	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1244	! -----------------------------------------------------------------
1245
1246	! Compute redistribution (advective) term
1247
1248	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1249	(rrd1omg(i,k)sigmaw(i)*d_th1(i,k) - &
1250	rrd2omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1)- &
1251	(1.-alpha_up(i,k))*omgbdth(i,k)- &
1252	alpha_up(i,k+1)omgbdth(i,k+1))ppi(i, k)
1253	! print*,'d_deltatw=', k, d_deltatw(i,k)
1254
1255	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1256	(rrd1omg(i,k)sigmaw(i)*d_q1(i,k)- &
1257	rrd2omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1)- &
1258	(1.-alpha_up(i,k))*omgbdq(i,k)- &
1259	alpha_up(i,k+1)*omgbdq(i,k+1))
1260	! print*,'d_deltaqw=', k, d_deltaqw(i,k)
1261
1262	! and increment large scale tendencies
1263
1264
1265
1266
1267	! C
1268	! -----------------------------------------------------------------
1269	d_tenv(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)- &
1270	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/ &
1271	(ph(i,k)-ph(i,k+1)) &
1272	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
1273	(ph(i,k)-ph(i,k+1)) )*ppi(i, k)
1274
1275	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)- &
1276	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/ &
1277	(ph(i,k)-ph(i,k+1)) &
1278	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
1279	(ph(i,k)-ph(i,k+1)) )
1280	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1281	d_tenv(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
1282
1283	d_qe(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
1284
1285	END IF
1286	! cc
1287	END DO
1288	END DO
1289	! ------------------------------------------------------------------
1290
1291	IF (prt_level>=10) THEN
1292	PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
1293	PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
1294	ENDIF
1295
1296	! Increment state variables
1297	!jyg<
1298	IF (iflag_wk_pop_dyn >= 1) THEN
1299	DO k = 1, klev
1300	DO i = 1, klon
1301	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1302	detr(i,k) = - d_sig_death(i) - d_sig_col(i)
1303	entr(i,k) = d_sig_gen(i)
1304	ENDIF
1305	ENDDO
1306	ENDDO
1307	ELSE ! (iflag_wk_pop_dyn >= 1)
1308	DO k = 1, klev
1309	DO i = 1, klon
1310	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1311	detr(i, k) = 0.
1312
1313	entr(i, k) = 0.
1314	ENDIF
1315	ENDDO
1316	ENDDO
1317	ENDIF ! (iflag_wk_pop_dyn >= 1)
1318
1319
1320
1321	DO k = 1, klev
1322	DO i = 1, klon
1323	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1324	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1325	! cc
1326
1327
1328
1329	! Coefficient de repartition
1330
1331	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1332	(ph(i,kupper(i))-ph(i,1))
1333	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
1334	(p(i,1)-ph(i,kupper(i)))
1335
1336
1337	! Reintroduce compensating subsidence term.
1338
1339	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1340	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1341	! . /(1-sigmaw)
1342	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1343	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1344	! . /(1-sigmaw)
1345
1346	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
1347	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
1348	! . /(1-sigmaw)
1349	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
1350	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
1351	! . /(1-sigmaw)
1352
1353	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1354	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1355	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1356
1357	!
1358
1359	! cc nrlmd Prise en compte du taux de mortalite
1360	! cc Definitions de entr, detr
1361	!jyg<
1362	!! detr(i, k) = 0.
1363	!!
1364	!! entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1365	!! sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1366	!!
1367	entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
1368	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1369	!>jyg
1370	wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1371
1372	! cc wkspread(i,k) =
1373	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
1374	! cc $ sigmaw(i)
1375
1376
1377	! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
1378	! Jingmei
1379
1380	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1381	! & Tgw(i,k),deltatw(i,k)
1382	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
1383	dtimesub
1384	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1385	ff(i) = d_deltatw(i, k)/dtimesub
1386
1387	! Sans GW
1388
1389	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-wkspread(k)deltatw(k))
1390
1391	! GW formule 1
1392
1393	! deltatw(k) = deltatw(k)+dtimesub*
1394	! $ (ff+dtKE(k) - wkspread(k)deltatw(k)-Tgw(k)deltatw(k))
1395
1396	! GW formule 2
1397
1398	IF (dtimesub*tgw(i,k)<1.E-10) THEN
1399	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
1400	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1401	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
1402	tgw(i,k)*deltatw(i,k) )
1403	ELSE
1404	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i,k)))* &
1405	(ff(i)+dtke(i,k) - &
1406	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1407	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - &
1408	tgw(i,k)*deltatw(i,k) )
1409	END IF
1410
1411	dth(i, k) = deltatw(i, k)/ppi(i, k)
1412
1413	gg(i) = d_deltaqw(i, k)/dtimesub
1414
1415	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
1416	entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
1417	(death_rate(i)sigmaw(i)+detr(i,k))deltaqw(i,k)/(1.-sigmaw(i)))
1418	! cc
1419
1420	! cc nrlmd
1421	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1422	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1423	! cc
1424	END IF
1425	END DO
1426	END DO
1427
1428
1429	! Scale tendencies so that water vapour remains positive in w and x.
1430
1431	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
1432	d_deltaqw, sigmaw, d_sigmaw, alpha)
1433	!
1434	! Alpha_tot = Product of all the alpha's
1435	DO i = 1, klon
1436	IF (wk_adv(i)) THEN
1437	alpha_tot(i) = alpha_tot(i)*alpha(i)
1438	END IF
1439	END DO
1440
1441	! cc nrlmd
1442	! c print*,'alpha'
1443	! c do i=1,klon
1444	! c print*,alpha(i)
1445	! c end do
1446	! cc
1447	DO k = 1, klev
1448	DO i = 1, klon
1449	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1450	d_tenv(i, k) = alpha(i)*d_tenv(i, k)
1451	d_qe(i, k) = alpha(i)*d_qe(i, k)
1452	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1453	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1454	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1455	END IF
1456	END DO
1457	END DO
1458	DO i = 1, klon
1459	IF (wk_adv(i)) THEN
1460	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1461	END IF
1462	END DO
1463
1464	! Update large scale variables and wake variables
1465	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1466	! IM 060208 DO k = 1,kupper(i)
1467	DO k = 1, klev
1468	DO i = 1, klon
1469	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1470	dtls(i, k) = dtls(i, k) + d_tenv(i, k)
1471	dqls(i, k) = dqls(i, k) + d_qe(i, k)
1472	! cc nrlmd
1473	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1474	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1475	! cc
1476	END IF
1477	END DO
1478	END DO
1479	DO k = 1, klev
1480	DO i = 1, klon
1481	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1482	tenv(i, k) = tenv0(i, k) + dtls(i, k)
1483	qe(i, k) = qe0(i, k) + dqls(i, k)
1484	the(i, k) = tenv(i, k)/ppi(i, k)
1485	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1486	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1487	dth(i, k) = deltatw(i, k)/ppi(i, k)
1488	! c print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
1489	! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1490	END IF
1491	END DO
1492	END DO
1493	!
1494	DO i = 1, klon
1495	IF (wk_adv(i)) THEN
1496	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1497	d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
1498	print *,'XXXX4 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
1499	END IF
1500	END DO
1501	!jyg<
1502	IF (iflag_wk_pop_dyn >= 1) THEN
1503	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1504	! Cumulatives
1505	DO i = 1, klon
1506	IF (wk_adv(i)) THEN
1507	d_sig_gen2(i) = d_sig_gen2(i) + d_sig_gen(i)
1508	d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
1509	d_sig_col2(i) = d_sig_col2(i) + d_sig_col(i)
1510	d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
1511	d_sig_bnd2(i) = d_sig_bnd2(i) + d_sig_bnd(i)
1512	END IF
1513	END DO
1514	! Bounds
1515	DO i = 1, klon
1516	IF (wk_adv(i)) THEN
1517	sigmaw_targ = max(sigmaw(i),sigmad)
1518	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
1519	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1520	print *,'XXXX5 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
1521	sigmaw(i) = sigmaw_targ
1522	END IF
1523	END DO
1524	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1525	! Cumulatives
1526	DO i = 1, klon
1527	IF (wk_adv(i)) THEN
1528	wdens(i) = wdens(i) + d_wdens(i)
1529	d_wdens2(i) = d_wdens2(i) + d_wdens(i)
1530	d_dens_gen2(i) = d_dens_gen2(i) + d_dens_gen(i)
1531	d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
1532	d_dens_col2(i) = d_dens_col2(i) + d_dens_col(i)
1533	d_dens_bnd2(i) = d_dens_bnd2(i) + d_dens_bnd(i)
1534	END IF
1535	END DO
1536	! Bounds
1537	DO i = 1, klon
1538	IF (wk_adv(i)) THEN
1539	wdens_targ = max(wdens(i),wdensmin)
1540	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
1541	d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
1542	wdens(i) = wdens_targ
1543	END IF
1544	END DO
1545	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1546	! Cumulatives
1547	DO i = 1, klon
1548	IF (wk_adv(i)) THEN
1549	awdens(i) = awdens(i) + d_awdens(i)
1550	d_awdens2(i) = d_awdens2(i) + d_awdens(i)
1551	END IF
1552	END DO
1553	! Bounds
1554	DO i = 1, klon
1555	IF (wk_adv(i)) THEN
1556	wdens_targ = min( max(awdens(i),0.), wdens(i) )
1557	d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
1558	awdens(i) = wdens_targ
1559	END IF
1560	END DO
1561	!
1562	IF (iflag_wk_pop_dyn == 2) THEN
1563	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn = 2!!!!!!
1564	! Cumulatives
1565	DO i = 1, klon
1566	IF (wk_adv(i)) THEN
1567	d_adens_death2(i) = d_adens_death2(i) + d_adens_death(i)
1568	d_adens_icol2(i) = d_adens_icol2(i) + d_adens_icol(i)
1569	d_adens_acol2(i) = d_adens_acol2(i) + d_adens_acol(i)
1570	d_adens_bnd2(i) = d_adens_bnd2(i) + d_adens_bnd(i)
1571	END IF
1572	END DO
1573	! Bounds
1574	DO i = 1, klon
1575	IF (wk_adv(i)) THEN
1576	wdens_targ = min( max(awdens(i),0.), wdens(i) )
1577	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
1578	END IF
1579	END DO
1580	ENDIF ! (iflag_wk_pop_dyn == 2)
1581	ENDIF ! (iflag_wk_pop_dyn >= 1)
1582
1583
1584	! Determine Ptop from buoyancy integral
1585	! ---------------------------------------
1586
1587	! - 1/ Pressure of the level where dth changes sign.
1588
1589	DO i = 1, klon
1590	IF (wk_adv(i)) THEN
1591	ptop_provis(i) = ph(i, 1)
1592	END IF
1593	END DO
1594
1595	DO k = 2, klev
1596	DO i = 1, klon
1597	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1598	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1599	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1600	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1601	END IF
1602	END DO
1603	END DO
1604
1605	! - 2/ dth integral
1606
1607	DO i = 1, klon
1608	IF (wk_adv(i)) THEN !!! nrlmd
1609	sum_dth(i) = 0.
1610	dthmin(i) = -delta_t_min
1611	z(i) = 0.
1612	END IF
1613	END DO
1614
1615	DO k = 1, klev
1616	DO i = 1, klon
1617	IF (wk_adv(i)) THEN
1618	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
1619	IF (dz(i)>0) THEN
1620	z(i) = z(i) + dz(i)
1621	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1622	dthmin(i) = amin1(dthmin(i), dth(i,k))
1623	END IF
1624	END IF
1625	END DO
1626	END DO
1627
1628	! - 3/ height of triangle with area= sum_dth and base = dthmin
1629
1630	DO i = 1, klon
1631	IF (wk_adv(i)) THEN
1632	hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1633	hw(i) = amax1(hwmin, hw(i))
1634	END IF
1635	END DO
1636
1637	! - 4/ now, get Ptop
1638
1639	DO i = 1, klon
1640	IF (wk_adv(i)) THEN !!! nrlmd
1641	ktop(i) = 0
1642	z(i) = 0.
1643	END IF
1644	END DO
1645
1646	DO k = 1, klev
1647	DO i = 1, klon
1648	IF (wk_adv(i)) THEN
1649	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw(i)-z(i))
1650	IF (dz(i)>0) THEN
1651	z(i) = z(i) + dz(i)
1652	ptop(i) = ph(i, k) - rho(i, k)RGdz(i)
1653	ktop(i) = k
1654	END IF
1655	END IF
1656	END DO
1657	END DO
1658
1659	! 4.5/Correct ktop and ptop
1660
1661	DO i = 1, klon
1662	IF (wk_adv(i)) THEN
1663	ptop_new(i) = ptop(i)
1664	END IF
1665	END DO
1666
1667	DO k = klev, 2, -1
1668	DO i = 1, klon
1669	! IM v3JYG; IF (k .GE. ktop(i)
1670	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1671	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1672	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1673	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1674	END IF
1675	END DO
1676	END DO
1677
1678
1679	DO i = 1, klon
1680	IF (wk_adv(i)) THEN
1681	ptop(i) = ptop_new(i)
1682	END IF
1683	END DO
1684
1685	DO k = klev, 1, -1
1686	DO i = 1, klon
1687	IF (wk_adv(i)) THEN !!! nrlmd
1688	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
1689	END IF
1690	END DO
1691	END DO
1692
1693	! 5/ Set deltatw & deltaqw to 0 above kupper
1694
1695	DO k = 1, klev
1696	DO i = 1, klon
1697	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1698	deltatw(i, k) = 0.
1699	deltaqw(i, k) = 0.
1700	d_deltatw2(i,k) = -deltatw0(i,k)
1701	d_deltaqw2(i,k) = -deltaqw0(i,k)
1702	END IF
1703	END DO
1704	END DO
1705
1706
1707	! -------------Cstar computation---------------------------------
1708	DO i = 1, klon
1709	IF (wk_adv(i)) THEN !!! nrlmd
1710	sum_thu(i) = 0.
1711	sum_tu(i) = 0.
1712	sum_qu(i) = 0.
1713	sum_thvu(i) = 0.
1714	sum_dth(i) = 0.
1715	sum_dq(i) = 0.
1716	sum_rho(i) = 0.
1717	sum_dtdwn(i) = 0.
1718	sum_dqdwn(i) = 0.
1719
1720	av_thu(i) = 0.
1721	av_tu(i) = 0.
1722	av_qu(i) = 0.
1723	av_thvu(i) = 0.
1724	av_dth(i) = 0.
1725	av_dq(i) = 0.
1726	av_rho(i) = 0.
1727	av_dtdwn(i) = 0.
1728	av_dqdwn(i) = 0.
1729	END IF
1730	END DO
1731
1732	! Integrals (and wake top level number)
1733	! --------------------------------------
1734
1735	! Initialize sum_thvu to 1st level virt. pot. temp.
1736
1737	DO i = 1, klon
1738	IF (wk_adv(i)) THEN !!! nrlmd
1739	z(i) = 1.
1740	dz(i) = 1.
1741	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
1742	sum_dth(i) = 0.
1743	END IF
1744	END DO
1745
1746	DO k = 1, klev
1747	DO i = 1, klon
1748	IF (wk_adv(i)) THEN !!! nrlmd
1749	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
1750	IF (dz(i)>0) THEN
1751	z(i) = z(i) + dz(i)
1752	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1753	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1754	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1755	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
1756	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1757	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1758	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1759	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1760	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1761	END IF
1762	END IF
1763	END DO
1764	END DO
1765
1766	DO i = 1, klon
1767	IF (wk_adv(i)) THEN !!! nrlmd
1768	hw0(i) = z(i)
1769	END IF
1770	END DO
1771
1772
1773	! - WAPE and mean forcing computation
1774	! ---------------------------------------
1775
1776	! ---------------------------------------
1777
1778	! Means
1779
1780	DO i = 1, klon
1781	IF (wk_adv(i)) THEN !!! nrlmd
1782	av_thu(i) = sum_thu(i)/hw0(i)
1783	av_tu(i) = sum_tu(i)/hw0(i)
1784	av_qu(i) = sum_qu(i)/hw0(i)
1785	av_thvu(i) = sum_thvu(i)/hw0(i)
1786	av_dth(i) = sum_dth(i)/hw0(i)
1787	av_dq(i) = sum_dq(i)/hw0(i)
1788	av_rho(i) = sum_rho(i)/hw0(i)
1789	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1790	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1791
1792	wape(i) = -RGhw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
1793	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
1794	END IF
1795	END DO
1796
1797	! Filter out bad wakes
1798
1799	DO k = 1, klev
1800	DO i = 1, klon
1801	IF (wk_adv(i)) THEN !!! nrlmd
1802	IF (wape(i)<0.) THEN
1803	deltatw(i, k) = 0.
1804	deltaqw(i, k) = 0.
1805	dth(i, k) = 0.
1806	d_deltatw2(i,k) = -deltatw0(i,k)
1807	d_deltaqw2(i,k) = -deltaqw0(i,k)
1808	END IF
1809	END IF
1810	END DO
1811	END DO
1812
1813	DO i = 1, klon
1814	IF (wk_adv(i)) THEN !!! nrlmd
1815	IF (wape(i)<0.) THEN
1816	wape(i) = 0.
1817	cstar(i) = 0.
1818	hw(i) = hwmin
1819	!jyg<
1820	!! sigmaw(i) = max(sigmad, sigd_con(i))
1821	sigmaw_targ = max(sigmad, sigd_con(i))
1822	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
1823	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1824	print *,'XXXX6 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
1825	sigmaw(i) = sigmaw_targ
1826	!>jyg
1827	fip(i) = 0.
1828	gwake(i) = .FALSE.
1829	ELSE
1830	cstar(i) = starksqrt(2.wape(i))
1831	gwake(i) = .TRUE.
1832	END IF
1833	END IF
1834	END DO
1835
1836	END DO ! end sub-timestep loop
1837
1838	IF (prt_level>=10) THEN
1839	PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
1840	sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
1841	ENDIF
1842
1843
1844	! ----------------------------------------------------------
1845	! Determine wake final state; recompute wape, cstar, ktop;
1846	! filter out bad wakes.
1847	! ----------------------------------------------------------
1848
1849	! 2.1 - Undisturbed area and Wake integrals
1850	! ---------------------------------------------------------
1851
1852	DO i = 1, klon
1853	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
1854	IF (ok_qx_qw(i)) THEN
1855	! cc
1856	z(i) = 0.
1857	sum_thu(i) = 0.
1858	sum_tu(i) = 0.
1859	sum_qu(i) = 0.
1860	sum_thvu(i) = 0.
1861	sum_dth(i) = 0.
1862	sum_half_dth(i) = 0.
1863	sum_dq(i) = 0.
1864	sum_rho(i) = 0.
1865	sum_dtdwn(i) = 0.
1866	sum_dqdwn(i) = 0.
1867
1868	av_thu(i) = 0.
1869	av_tu(i) = 0.
1870	av_qu(i) = 0.
1871	av_thvu(i) = 0.
1872	av_dth(i) = 0.
1873	av_dq(i) = 0.
1874	av_rho(i) = 0.
1875	av_dtdwn(i) = 0.
1876	av_dqdwn(i) = 0.
1877
1878	dthmin(i) = -delta_t_min
1879	END IF
1880	END DO
1881	! Potential temperatures and humidity
1882	! ----------------------------------------------------------
1883
1884	DO k = 1, klev
1885	DO i = 1, klon
1886	! cc nrlmd IF ( wk_adv(i)) THEN
1887	IF (ok_qx_qw(i)) THEN
1888	! cc
1889	rho(i, k) = p(i, k)/(RD*tenv(i,k))
1890	IF (k==1) THEN
1891	rhoh(i, k) = ph(i, k)/(RD*tenv(i,k))
1892	zhh(i, k) = 0
1893	ELSE
1894	rhoh(i, k) = ph(i, k)2./(RD(tenv(i,k)+tenv(i,k-1)))
1895	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
1896	END IF
1897	the(i, k) = tenv(i, k)/ppi(i, k)
1898	thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1899	tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i)
1900	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
1901	rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k)))
1902	dth(i, k) = deltatw(i, k)/ppi(i, k)
1903	END IF
1904	END DO
1905	END DO
1906
1907	! Integrals (and wake top level number)
1908	! -----------------------------------------------------------
1909
1910	! Initialize sum_thvu to 1st level virt. pot. temp.
1911
1912	DO i = 1, klon
1913	! cc nrlmd IF ( wk_adv(i)) THEN
1914	IF (ok_qx_qw(i)) THEN
1915	! cc
1916	z(i) = 1.
1917	dz(i) = 1.
1918	dz_half(i) = 1.
1919	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
1920	sum_dth(i) = 0.
1921	END IF
1922	END DO
1923
1924	DO k = 1, klev
1925	DO i = 1, klon
1926	! cc nrlmd IF ( wk_adv(i)) THEN
1927	IF (ok_qx_qw(i)) THEN
1928	! cc
1929	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
1930	dz_half(i) = -(amax1(ph(i,k+1),0.5(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)RG)
1931	IF (dz(i)>0) THEN
1932	z(i) = z(i) + dz(i)
1933	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1934	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1935	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1936	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
1937	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1938	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1939	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1940	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1941	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1942	!
1943	dthmin(i) = min(dthmin(i), dth(i,k))
1944	END IF
1945	IF (dz_half(i)>0) THEN
1946	sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
1947	END IF
1948	END IF
1949	END DO
1950	END DO
1951
1952	DO i = 1, klon
1953	! cc nrlmd IF ( wk_adv(i)) THEN
1954	IF (ok_qx_qw(i)) THEN
1955	! cc
1956	hw0(i) = z(i)
1957	END IF
1958	END DO
1959
1960	! - WAPE and mean forcing computation
1961	! -------------------------------------------------------------
1962
1963	! Means
1964
1965	DO i = 1, klon
1966	! cc nrlmd IF ( wk_adv(i)) THEN
1967	IF (ok_qx_qw(i)) THEN
1968	! cc
1969	av_thu(i) = sum_thu(i)/hw0(i)
1970	av_tu(i) = sum_tu(i)/hw0(i)
1971	av_qu(i) = sum_qu(i)/hw0(i)
1972	av_thvu(i) = sum_thvu(i)/hw0(i)
1973	av_dth(i) = sum_dth(i)/hw0(i)
1974	av_dq(i) = sum_dq(i)/hw0(i)
1975	av_rho(i) = sum_rho(i)/hw0(i)
1976	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1977	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1978
1979	wape2(i) = -RGhw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
1980	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
1981	END IF
1982	END DO
1983
1984
1985
1986	! Prognostic variable update
1987	! ------------------------------------------------------------
1988
1989	! Filter out bad wakes
1990
1991	IF (iflag_wk_check_trgl>=1) THEN
1992	! Check triangular shape of dth profile
1993	DO i = 1, klon
1994	IF (ok_qx_qw(i)) THEN
1995	!! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
1996	!! print ,'wake, 2.sum_dth(i)/(hw0(i)*dthmin(i)) ', &
1997	!! 2.sum_dth(i)/(hw0(i)dthmin(i))
1998	!! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
1999	!! sum_half_dth(i), sum_dth(i)
2000	IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
2001	wape2(i) = -1.
2002	!! print *,'wake, rej 1'
2003	ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.sum_dth(i)/(hw0(i)dthmin(i)) - 1.) > 0.5) THEN
2004	wape2(i) = -1.
2005	!! print *,'wake, rej 2'
2006	ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
2007	wape2(i) = -1.
2008	!! print *,'wake, rej 3'
2009	END IF
2010	END IF
2011	END DO
2012	END IF
2013
2014
2015	DO k = 1, klev
2016	DO i = 1, klon
2017	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
2018	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
2019	! cc
2020	deltatw(i, k) = 0.
2021	deltaqw(i, k) = 0.
2022	dth(i, k) = 0.
2023	d_deltatw2(i,k) = -deltatw0(i,k)
2024	d_deltaqw2(i,k) = -deltaqw0(i,k)
2025	END IF
2026	END DO
2027	END DO
2028
2029
2030	DO i = 1, klon
2031	! cc nrlmd IF ( wk_adv(i)) THEN
2032	IF (ok_qx_qw(i)) THEN
2033	! cc
2034	IF (wape2(i)<0.) THEN
2035	wape2(i) = 0.
2036	cstar2(i) = 0.
2037	hw(i) = hwmin
2038	!jyg<
2039	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
2040	sigmaw_targ = max(sigmad, sigd_con(i))
2041	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
2042	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2043	print *,'XXXX7 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
2044	sigmaw(i) = sigmaw_targ
2045	!>jyg
2046	fip(i) = 0.
2047	gwake(i) = .FALSE.
2048	ELSE
2049	IF (prt_level>=10) PRINT *, 'wape2>0'
2050	cstar2(i) = starksqrt(2.wape2(i))
2051	gwake(i) = .TRUE.
2052	END IF
2053	END IF
2054	END DO
2055
2056	DO i = 1, klon
2057	! cc nrlmd IF ( wk_adv(i)) THEN
2058	IF (ok_qx_qw(i)) THEN
2059	! cc
2060	ktopw(i) = ktop(i)
2061	END IF
2062	END DO
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 (ktopw(i)>0 .AND. gwake(i)) THEN
2069
2070	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
2071	! cc heff = 600.
2072	! Utilisation de la hauteur hw
2073	! c heff = 0.7*hw
2074	heff(i) = hw(i)
2075
2076	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
2077	sqrt(sigmaw(i)wdens(i)3.14)
2078	fip(i) = alpk*fip(i)
2079	! jyg2
2080	ELSE
2081	fip(i) = 0.
2082	END IF
2083	END IF
2084	END DO
2085
2086	! Limitation de sigmaw
2087
2088	! cc nrlmd
2089	! DO i=1,klon
2090	! IF (OK_qx_qw(i)) THEN
2091	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
2092	! ENDIF
2093	! ENDDO
2094	! cc
2095
2096	!jyg<
2097	IF (iflag_wk_pop_dyn >= 1) THEN
2098	DO i = 1, klon
2099	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2100	.NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
2101	ENDDO
2102	ELSE ! (iflag_wk_pop_dyn >= 1)
2103	DO i = 1, klon
2104	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2105	.NOT. ok_qx_qw(i)
2106	ENDDO
2107	ENDIF ! (iflag_wk_pop_dyn >= 1)
2108	!>jyg
2109
2110	DO k = 1, klev
2111	DO i = 1, klon
2112	!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2113	!!jyg .NOT. ok_qx_qw(i)) THEN
2114	IF (kill_wake(i)) THEN
2115	! cc
2116	dtls(i, k) = 0.
2117	dqls(i, k) = 0.
2118	deltatw(i, k) = 0.
2119	deltaqw(i, k) = 0.
2120	d_deltatw2(i,k) = -deltatw0(i,k)
2121	d_deltaqw2(i,k) = -deltaqw0(i,k)
2122	END IF ! (kill_wake(i))
2123	END DO
2124	END DO
2125
2126	DO i = 1, klon
2127	!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2128	!!jyg .NOT. ok_qx_qw(i)) THEN
2129	IF (kill_wake(i)) THEN
2130	ktopw(i) = 0
2131	wape(i) = 0.
2132	cstar(i) = 0.
2133	!!jyg Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
2134	!! hw(i) = hwmin !jyg
2135	!! sigmaw(i) = sigmad !jyg
2136	hw(i) = 0. !jyg
2137	fip(i) = 0.
2138	!! sigmaw(i) = 0. !jyg
2139	sigmaw_targ = 0.
2140	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
2141	!! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2142	d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i) ! _in = correction jyg 20220124
2143	print *,'XXXX8 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)
2144	sigmaw(i) = sigmaw_targ
2145	IF (iflag_wk_pop_dyn >= 1) THEN
2146	!! awdens(i) = 0.
2147	!! wdens(i) = 0.
2148	wdens_targ = 0.
2149	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
2150	!! d_wdens2(i) = wdens_targ - wdens(i)
2151	d_wdens2(i) = wdens_targ - wdens_in(i) ! jyg 20220916
2152	wdens(i) = wdens_targ
2153	wdens_targ = 0.
2154	!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
2155	IF (iflag_wk_pop_dyn == 2) THEN
2156	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
2157	ENDIF ! (iflag_wk_pop_dyn == 2)
2158	!! d_awdens2(i) = wdens_targ - awdens(i)
2159	d_awdens2(i) = wdens_targ - awdens_in(i) ! jyg 20220916
2160	awdens(i) = wdens_targ
2161	!! IF (iflag_wk_pop_dyn == 2) THEN
2162	!! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
2163	!! ENDIF ! (iflag_wk_pop_dyn == 2)
2164	ENDIF ! (iflag_wk_pop_dyn >= 1)
2165	ELSE ! (kill_wake(i))
2166	wape(i) = wape2(i)
2167	cstar(i) = cstar2(i)
2168	END IF ! (kill_wake(i))
2169	! c print*,'wape wape2 ktopw OK_qx_qw =',
2170	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
2171	END DO
2172
2173	IF (prt_level>=10) THEN
2174	PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
2175	wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
2176	ENDIF
2177
2178
2179	! -----------------------------------------------------------------
2180	! Get back to tendencies per second
2181
2182	DO k = 1, klev
2183	DO i = 1, klon
2184
2185	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
2186	!jyg<
2187	!! IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
2188	IF (ok_qx_qw(i)) THEN
2189	!>jyg
2190	! cc
2191	dtls(i, k) = dtls(i, k)/dtime
2192	dqls(i, k) = dqls(i, k)/dtime
2193	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
2194	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
2195	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
2196	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
2197	! c $ ,death_rate(i)*sigmaw(i)
2198	END IF
2199	END DO
2200	END DO
2201	!jyg<
2202	IF (iflag_wk_pop_dyn >= 1) THEN
2203	DO i = 1, klon
2204	IF (ok_qx_qw(i)) THEN
2205	d_sig_gen2(i) = d_sig_gen2(i)/dtime
2206	d_sig_death2(i) = d_sig_death2(i)/dtime
2207	d_sig_col2(i) = d_sig_col2(i)/dtime
2208	d_sig_spread2(i) = d_sig_spread2(i)/dtime
2209	d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
2210	d_sigmaw2(i) = d_sigmaw2(i)/dtime
2211	print *,'XXXX9 d_sigmaw2(i), sigmaw(i), dtime ', d_sigmaw2(i), sigmaw(i), dtime
2212	!
2213	d_dens_gen2(i) = d_dens_gen2(i)/dtime
2214	d_dens_death2(i) = d_dens_death2(i)/dtime
2215	d_dens_col2(i) = d_dens_col2(i)/dtime
2216	d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
2217	d_awdens2(i) = d_awdens2(i)/dtime
2218	d_wdens2(i) = d_wdens2(i)/dtime
2219	ENDIF
2220	ENDDO
2221	IF (iflag_wk_pop_dyn == 2) THEN
2222	DO i = 1, klon
2223	IF (ok_qx_qw(i)) THEN
2224	d_adens_death2(i) = d_adens_death2(i)/dtime
2225	d_adens_icol2(i) = d_adens_icol2(i)/dtime
2226	d_adens_acol2(i) = d_adens_acol2(i)/dtime
2227	d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
2228	ENDIF
2229	ENDDO
2230	ENDIF ! (iflag_wk_pop_dyn == 2)
2231	ENDIF ! (iflag_wk_pop_dyn >= 1)
2232
2233	!>jyg
2234
2235	RETURN
2236	END SUBROUTINE wake
2237
2238	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
2239	d_deltaqw, sigmaw, d_sigmaw, alpha)
2240	! ------------------------------------------------------
2241	! Dtermination du coefficient alpha tel que les tendances
2242	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
2243	! a une humidite positive dans la zone (x) et dans la zone (w).
2244	! ------------------------------------------------------
2245	IMPLICIT NONE
2246
2247	! Input
2248	REAL qe(nlon, nl), d_qe(nlon, nl)
2249	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
2250	REAL sigmaw(nlon), d_sigmaw(nlon)
2251	LOGICAL wk_adv(nlon)
2252	INTEGER nl, nlon
2253	! Output
2254	REAL alpha(nlon)
2255	! Internal variables
2256	REAL zeta(nlon, nl)
2257	REAL alpha1(nlon)
2258	REAL x, a, b, c, discrim
2259	REAL epsilon_loc
2260	INTEGER i,k
2261
2262	DO k = 1, nl
2263	DO i = 1, nlon
2264	IF (wk_adv(i)) THEN
2265	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
2266	zeta(i, k) = 0.
2267	ELSE
2268	zeta(i, k) = 1.
2269	END IF
2270	END IF
2271	END DO
2272	DO i = 1, nlon
2273	IF (wk_adv(i)) THEN
2274	x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
2275	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i) &
2276	(deltaqw(i,k)+d_deltaqw(i,k))
2277	a = -d_sigmaw(i)*d_deltaqw(i, k)
2278	b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
2279	deltaqw(i, k)*d_sigmaw(i)
2280	c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc
2281	discrim = bb - 4.a*c
2282	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
2283	IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap
2284	alpha1(i) = 1.
2285	ELSE
2286	IF (x>=0.) THEN
2287	alpha1(i) = 1.
2288	ELSE
2289	IF (a>0.) THEN
2290	alpha1(i) = 0.9min( (2.c)/(-b+sqrt(discrim)), &
2291	(-b+sqrt(discrim))/(2.*a) )
2292	ELSE IF (a==0.) THEN
2293	alpha1(i) = 0.9*(-c/b)
2294	ELSE
2295	! print*,'a,b,c discrim',a,b,c discrim
2296	alpha1(i) = 0.9max( (2.c)/(-b+sqrt(discrim)), &
2297	(-b+sqrt(discrim))/(2.*a))
2298	END IF
2299	END IF
2300	END IF
2301	alpha(i) = min(alpha(i), alpha1(i))
2302	END IF
2303	END DO
2304	END DO
2305
2306	RETURN
2307	END SUBROUTINE wake_vec_modulation
2308
2309
2310
2311	SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper)
2312
2313	USE wake_ini_mod , ONLY : pupperbyphs
2314	IMPLICIT NONE
2315
2316	INTEGER, INTENT(IN) :: klon,klev
2317	REAL, INTENT(IN), DIMENSION (klon,klev+1) :: ph
2318	REAL, INTENT(IN), DIMENSION (klon) :: ptop
2319	REAL, INTENT(OUT), DIMENSION (klon) :: pupper
2320	INTEGER, INTENT(OUT), DIMENSION (klon) :: kupper
2321	INTEGER :: i,k
2322
2323
2324	kupper = 0
2325
2326	IF (pupperbyphs<1.) THEN
2327	! Choose an integration bound well above wake top
2328	! -----------------------------------------------------------------
2329
2330	! Pupper = 50000. ! melting level
2331	! Pupper = 60000.
2332	! Pupper = 80000. ! essais pour case_e
2333	DO i = 1, klon
2334	! pupper(i) = 0.6*ph(i, 1)
2335	pupper(i) = pupperbyphs*ph(i, 1)
2336	pupper(i) = max(pupper(i), 45000.)
2337	! cc Pupper(i) = 60000.
2338	END DO
2339
2340	ELSE
2341
2342	DO i=1, klon
2343	! pupper(i) = pupperbyphsptop(i)+(1.-pupperbyphs)ph(i, 1)
2344	pupper(i) = min( pupperbyphsptop(i)+(1.-pupperbyphs)ph(i, 1) , ptop(i)-5000.)
2345	END DO
2346	END IF
2347
2348	! -5/ Determination de kupper
2349
2350	DO k = klev, 1, -1
2351	DO i = 1, klon
2352	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
2353	END DO
2354	END DO
2355
2356	! On evite kupper = 1 et kupper = klev
2357	DO i = 1, klon
2358	kupper(i) = max(kupper(i), 2)
2359	kupper(i) = min(kupper(i), klev-1)
2360	END DO
2361	RETURN
2362	END SUBROUTINE pkupper
2363
2364
2365	SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
2366	dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
2367	d_awdens, d_wdens, d_sigmaw, &
2368	iflag_wk_act, wk_adv, cin, wape, &
2369	drdt, &
2370	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
2371	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
2372	d_wdens_targ, d_sigmaw_targ)
2373
2374
2375	USE wake_ini_mod , ONLY : wake_ini
2376	USE wake_ini_mod , ONLY : prt_level,RG
2377	USE wake_ini_mod , ONLY : stark, wdens_ref
2378	USE wake_ini_mod , ONLY : tau_cv, rzero, aa0
2379	USE wake_ini_mod , ONLY : iflag_wk_pop_dyn, wdensmin
2380	USE wake_ini_mod , ONLY : sigmad, cstart, sigmaw_max
2381
2382	IMPLICIT NONE
2383
2384	INTEGER, INTENT(IN) :: klon,klev
2385	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
2386	REAL, INTENT(IN) :: dtime
2387	REAL, INTENT(IN) :: dtimesub
2388	REAL, DIMENSION (klon), INTENT(IN) :: wgen
2389	REAL, DIMENSION (klon), INTENT(IN) :: wdens
2390	REAL, DIMENSION (klon), INTENT(IN) :: awdens
2391	REAL, DIMENSION (klon), INTENT(IN) :: sigmaw
2392	REAL, DIMENSION (klon), INTENT(IN) :: gfl, cstar
2393	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape
2394	REAL, DIMENSION (klon), INTENT(IN) :: rad_wk
2395	REAL, DIMENSION (klon), INTENT(IN) :: f_shear
2396	INTEGER, INTENT(IN) :: iflag_wk_act
2397
2398
2399	!
2400
2401	! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
2402	! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
2403	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_awdens, d_wdens
2404	REAL, DIMENSION (klon), INTENT(OUT) :: drdt
2405	! Some components of the tendencies of state variables
2406	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
2407	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_spread
2408	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
2409	REAL, INTENT(OUT) :: d_wdens_targ, d_sigmaw_targ
2410
2411
2412	REAL :: delta_t_min
2413	INTEGER :: nsub
2414	INTEGER :: i, k
2415	REAL :: wdens0
2416	! IM 080208
2417	LOGICAL, DIMENSION (klon) :: gwake
2418
2419	! Variables liees a la dynamique de population
2420	REAL, DIMENSION(klon) :: act
2421	REAL, DIMENSION(klon) :: tau_wk_inv
2422	REAL, DIMENSION(klon) :: wape1_act, wape2_act
2423	LOGICAL, DIMENSION (klon) :: kill_wake
2424	REAL :: drdt_pos
2425	REAL :: tau_wk_inv_min
2426
2427
2428
2429	IF (iflag_wk_act == 0) THEN
2430	act(:) = 0.
2431	ELSEIF (iflag_wk_act == 1) THEN
2432	act(:) = 1.
2433	ELSEIF (iflag_wk_act ==2) THEN
2434	DO i = 1, klon
2435	IF (wk_adv(i)) THEN
2436	wape1_act(i) = abs(cin(i))
2437	wape2_act(i) = 2.*wape1_act(i) + 1.
2438	act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
2439	ENDIF ! (wk_adv(i))
2440	ENDDO
2441	ENDIF ! (iflag_wk_act ==2)
2442
2443
2444	DO i = 1, klon
2445	! print*, 'XXX wk_adv(i)', wk_adv(i)
2446	IF (wk_adv(i)) THEN
2447	!! tau_wk(i) = max(rad_wk(i)/(3.cstar(i))((cstar(i)/cstart)**1.5 - 1), 100.)
2448	tau_wk_inv(i) = max( (3.cstar(i))/(rad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
2449	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
2450	drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
2451	(1 + 2f_shear(i)(2.sigmaw(i)-aa0wdens(i)) - 2.*sigmaw(i))
2452	!! (1 - 2sigmaw(i)(1.-f_shear(i)))
2453	drdt_pos=max(drdt(i),0.)
2454
2455	!! d_wdens(i) = ( wgen(i)(1.+2.(sigmaw(i)-sigmad)) &
2456	!! - wdens(i)*tau_wk_inv_min &
2457	!! - 2.gfl(i)wdens(i)Cstar(i) )dtimesub
2458	!jyg+mlt<
2459	d_awdens(i) = ( wgen(i) - (1./tau_cv)(awdens(i) - act(i)wdens(i)) )*dtimesub
2460	d_dens_gen(i) = wgen(i)
2461	d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
2462	d_dens_col(i) = -2.wdens(i)gfl(i)*drdt_pos
2463	d_dens_gen(i) = d_dens_gen(i)*dtimesub
2464	d_dens_death(i) = d_dens_death(i)*dtimesub
2465	d_dens_col(i) = d_dens_col(i)*dtimesub
2466
2467	d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
2468	!! d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min - &
2469	!! 2.wdens(i)gfl(i)drdt_pos )dtimesub
2470	!>jyg+mlt
2471	!
2472	!jyg<
2473	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
2474	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
2475	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
2476	d_wdens(i) = d_wdens_targ
2477	!! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
2478	!>jyg
2479
2480	!jyg+mlt<
2481	!! d_sigmaw(i) = ( (1.-2f_shear(i)sigmaw(i))(gfl(i)Cstar(i)+wgen(i)*sigmad/wdens(i)) &
2482	!! + 2.f_shear(i)wgen(i)sigmaw(i)*2/wdens(i) &
2483	!! - sigmaw(i)tau_wk_inv_min )dtimesub
2484	d_sig_gen(i) = wgen(i)*aa0
2485	!! print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', &
2486	!! sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min
2487	d_sig_death(i) = - sigmaw(i)(1.-awdens(i)/wdens(i))tau_wk_inv_min
2488	!!
2489
2490	d_sig_col(i) = - 2f_shear(i)sigmaw(i)gfl(i)drdt_pos
2491	d_sig_col(i) = - 2f_shear(i)(2.sigmaw(i)-wdens(i)aa0)gfl(i)drdt_pos
2492	d_sig_spread(i) = gfl(i)*cstar(i)
2493	d_sig_gen(i) = d_sig_gen(i)*dtimesub
2494	d_sig_death(i) = d_sig_death(i)*dtimesub
2495	d_sig_col(i) = d_sig_col(i)*dtimesub
2496	d_sig_spread(i) = d_sig_spread(i)*dtimesub
2497	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
2498	!>jyg+mlt
2499	!
2500	!jyg<
2501	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
2502	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
2503	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
2504	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
2505	d_sigmaw(i) = d_sigmaw_targ
2506	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
2507	!>jyg
2508	ENDIF
2509	ENDDO
2510
2511	IF (prt_level >= 10) THEN
2512	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
2513	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
2514	print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
2515	wdens(1), awdens(1), act(1), d_awdens(1)
2516	print ,'wake, wgen, -(wdens-awdens)tau_wk_inv, -2.wdensgfl*drdt_pos, d_wdens ', &
2517	wgen(1), -(wdens(1)-awdens(1))tau_wk_inv(1), -2.wdens(1)gfl(1)drdt_pos, d_wdens(1)
2518	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
2519	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
2520	ENDIF
2521
2522	RETURN
2523	END SUBROUTINE wake_popdyn_1
2524
2525	SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
2526	sigmaw, wdens, awdens, & !! states variables
2527	gfl, cstar, cin, wape, rad_wk, &
2528	d_sigmaw, d_wdens, d_awdens, & !! tendences
2529	cont_fact, &
2530	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
2531	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
2532	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
2533
2534
2535
2536	USE wake_ini_mod , ONLY : wake_ini
2537	USE wake_ini_mod , ONLY : prt_level,RG
2538	USE wake_ini_mod , ONLY : stark, wdens_ref
2539	USE wake_ini_mod , ONLY : tau_cv, rzero, aa0
2540	USE wake_ini_mod , ONLY : iflag_wk_pop_dyn, wdensmin
2541	USE wake_ini_mod , ONLY : sigmad, cstart, sigmaw_max
2542
2543	IMPLICIT NONE
2544
2545	INTEGER, INTENT(IN) :: klon,klev
2546	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
2547	REAL, INTENT(IN) :: dtimesub
2548	REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes
2549	REAL, DIMENSION (klon), INTENT(IN) :: sigmaw !! sigma = fractional area of wakes
2550	REAL, DIMENSION (klon), INTENT(IN) :: wdens !! D = number of wakes per unit area
2551	REAL, DIMENSION (klon), INTENT(IN) :: awdens !! A = number of active wakes per unit area
2552	REAL, DIMENSION (klon), INTENT(IN) :: gfl !! Lg = gust front lenght per unit area
2553	REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes
2554	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape
2555	REAL, DIMENSION (klon), INTENT(IN) :: rad_wk !! r = wake radius
2556
2557	!
2558	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens
2559	REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact
2560	! Some components of the tendencies of state variables
2561	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
2562	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
2563	REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
2564
2565
2566	!! internal variables
2567
2568	INTEGER :: i, k
2569	REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes
2570	REAL :: tau_wk_inv_min
2571	REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes
2572	REAL :: d_wdens_targ, d_sigmaw_targ
2573
2574
2575	!! Equations
2576	!! dD/dt = B - (D-A)/tau - f D^2
2577	!! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2
2578	!! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0)
2579	!!
2580	!! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma))
2581
2582
2583	DO i = 1, klon
2584	! print*, 'XXX wk_adv(i)', wk_adv(i)
2585	IF (wk_adv(i)) THEN
2586	!! tau_wk(i) = max(rad_wk(i)/(3.cstar(i))((cstar(i)/cstart)**1.5 - 1), 100.)
2587	tau_wk_inv(i) = max( (3.cstar(i))/(rad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
2588	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
2589	tau_prime(i) = tau_cv
2590	!! cont_fact(i) = 2.(wgen(i)(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / &
2591	!! (2.(wdens(i)-awdens(i))2(2.sigmaw(i)/wdens(i) - aa0) + wdens(i)(1.-2.*sigmaw(i)))
2592	cont_fact(i) = 2.3.14rad_wk(i)*cstar(i)
2593
2594	d_sig_gen(i) = wgen(i)*aa0
2595	d_sig_death(i) = - sigmaw(i)(1.-awdens(i)/wdens(i))tau_wk_inv_min
2596	d_sig_col(i) = - cont_fact(i)(wdens(i)-awdens(i))2(2.*sigmaw(i)/wdens(i)-aa0)
2597	d_sig_spread(i) = gfl(i)*cstar(i)
2598	!
2599	d_sig_gen(i) = d_sig_gen(i)*dtimesub
2600	d_sig_death(i) = d_sig_death(i)*dtimesub
2601	d_sig_col(i) = d_sig_col(i)*dtimesub
2602	d_sig_spread(i) = d_sig_spread(i)*dtimesub
2603	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
2604
2605
2606	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
2607	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
2608	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
2609	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
2610	d_sigmaw(i) = d_sigmaw_targ
2611	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
2612
2613
2614	d_dens_gen(i) = wgen(i)
2615	d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
2616	d_dens_col(i) = - cont_fact(i)wdens(i)*2
2617	!
2618	d_dens_gen(i) = d_dens_gen(i)*dtimesub
2619	d_dens_death(i) = d_dens_death(i)*dtimesub
2620	d_dens_col(i) = d_dens_col(i)*dtimesub
2621	d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
2622
2623
2624	d_adens_death(i) = -awdens(i)/tau_prime(i)
2625	d_adens_icol(i) = cont_fact(i)(wdens(i)-awdens(i))*2
2626	d_adens_acol(i) = - cont_fact(i)awdens(i)*2
2627	!
2628	d_adens_death(i) = d_adens_death(i)*dtimesub
2629	d_adens_icol(i) = d_adens_icol(i)*dtimesub
2630	d_adens_acol(i) = d_adens_acol(i)*dtimesub
2631	d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)
2632
2633	!!
2634	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
2635	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
2636	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
2637	d_wdens(i) = d_wdens_targ
2638
2639	d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
2640	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
2641	d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
2642	d_awdens(i) = d_wdens_targ
2643
2644
2645
2646	ENDIF
2647	ENDDO
2648
2649	IF (prt_level >= 10) THEN
2650	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', &
2651	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1)
2652	print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
2653	wdens(1), awdens(1), d_awdens(1)
2654	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
2655	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
2656	ENDIF
2657
2658	RETURN
2659	END SUBROUTINE wake_popdyn_2
2660

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

Download in other formats:

Original Format