Context Navigation

wake.F90 @ 3605

Last change on this file since 3605 was 3605, checked in by lguez, 5 years ago
Merge revisions 3427:3600 of trunk into branch `Ocean_skin`
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: 71.7 KB

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

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

Download in other formats:

Original Format