Context Navigation

wake.F90 @ 5440

Last change on this file since 5440 was 2787, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2727:2785 into testing branch
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: 62.0 KB

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

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

Download in other formats:

Original Format