Context Navigation

wake.F90 @ 2679

Last change on this file since 2679 was 2671, checked in by jyg, 8 years ago
Bug fix in cv3_routine.F90 (cv3_unsat could create a precipitating downdraught when convection was off) and in wake.F90 and calwake.F90 (some array dimensions were incompatible and some initializations were missing)
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: 59.5 KB

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

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

Download in other formats:

Original Format