Context Navigation

lmdz_wake.F90 @ 5101

Last change on this file since 5101 was 5101, checked in by abarral, 2 months ago
Handle DEBUG_IO in lmdz_cppkeys_wrapper.F90 Transform some files .F -> .[fF]90 [ne compile pas à cause de writefield_u non défini - en attente de réponse Laurent]
File size: 82.1 KB

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

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

Download in other formats:

Original Format