Context Navigation

wake.F90 @ 2308

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

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

Download in other formats:

Original Format