Context Navigation

wake.F90 @ 5017

Last change on this file since 5017 was 2495, checked in by lguez, 9 years ago

Bug fix. Virtual temperature is:

[1 + (1 / eps - 1) q] T

not:

(1 + eps q) T

(no worry: the numerical difference is small).

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: 53.0 KB

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

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

Download in other formats:

Original Format