Context Navigation

wake.F90 @ 5408

Last change on this file since 5408 was 3648, checked in by jghattas, 5 years ago
Correction d'erreur introduit en commit [3641] pour faire tourner en mode debug.
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: 72.2 KB

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

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

Download in other formats:

Original Format