Context Navigation

wake.F90 @ 3641

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

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

Download in other formats:

Original Format