Context Navigation

wake.F90 @ 2408

Last change on this file since 2408 was 2408, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2298:2396 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 53.0 KB

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

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

Download in other formats:

Original Format