Context Navigation

wake.F90 @ 3983

Last change on this file since 3983 was 3831, checked in by ymipsl, 10 years ago
module reorganisation for a cleaner dyn-phys interface YM
File size: 52.1 KB

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

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

Download in other formats:

Original Format