Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

wake.F90 @ 409

Last change on this file since 409 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 90.2 KB

Line
1	!
2	! $Id: wake.F 1669 2012-10-16 12:41:50Z fairhead $
3	!
4	SUBROUTINE WAKE (p,ph,pi,dtime,sigd_con &
5	& ,te0,qe0,omgb &
6	& ,dtdwn,dqdwn,amdwn,amup,dta,dqa &
7	& ,wdtPBL,wdqPBL,udtPBL,udqPBL &
8	& ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl &
9	& ,dtls,dqls &
10	& ,ktopw,omgbdth,dp_omgb,wdens &
11	& ,tu,qu &
12	& ,dtKE,dqKE &
13	& ,dtPBL,dqPBL &
14	& ,omg,dp_deltomg,spread &
15	& ,Cstar,d_deltat_gw &
16	& ,d_deltatw2,d_deltaqw2)
17
18
19	!***************************************************************
20	!* *
21	!* WAKE *
22	!* retour a un Pupper fixe *
23	!* *
24	!* written by : GRANDPEIX Jean-Yves 09/03/2000 *
25	!* modified by : ROEHRIG Romain 01/29/2007 *
26	!***************************************************************
27	!c
28	use dimphy
29	IMPLICIT none
30	!c============================================================================
31	!C
32	!C
33	!C But : Decrire le comportement des poches froides apparaissant dans les
34	!C grands systemes convectifs, et fournir l'energie disponible pour
35	!C le declenchement de nouvelles colonnes convectives.
36	!C
37	!C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area
38	!C deltaqw : ecart d'humidite wake-undisturbed area
39	!C sigmaw : fraction d'aire occupee par la poche.
40	!C
41	!C Variable de sortie :
42	!c
43	!c wape : WAke Potential Energy
44	!c fip : Front Incident Power (W/m2) - ALP
45	!c gfl : Gust Front Length per unit area (m-1)
46	!C dtls : large scale temperature tendency due to wake
47	!C dqls : large scale humidity tendency due to wake
48	!C hw : hauteur de la poche
49	!C dp_omgb : vertical gradient of large scale omega
50	!C wdens : densite de poches
51	!C omgbdth: flux of Delta_Theta transported by LS omega
52	!C dtKE : differential heating (wake - unpertubed)
53	!C dqKE : differential moistening (wake - unpertubed)
54	!C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
55	!C dp_deltomg : vertical gradient of omg (s-1)
56	!C spread : spreading term in dt_wake and dq_wake
57	!C deltatw : updated temperature difference (T_w-T_u).
58	!C deltaqw : updated humidity difference (q_w-q_u).
59	!C sigmaw : updated wake fractional area.
60	!C d_deltat_gw : delta T tendency due to GW
61	!c
62	!C Variables d'entree :
63	!c
64	!c aire : aire de la maille
65	!c te0 : temperature dans l'environnement (K)
66	!C qe0 : humidite dans l'environnement (kg/kg)
67	!C omgb : vitesse verticale moyenne sur la maille (Pa/s)
68	!C dtdwn: source de chaleur due aux descentes (K/s)
69	!C dqdwn: source d'humidite due aux descentes (kg/kg/s)
70	!C dta : source de chaleur due courants satures et detrain (K/s)
71	!C dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
72	!C amdwn: flux de masse total des descentes, par unite de
73	!C surface de la maille (kg/m2/s)
74	!C amup : flux de masse total des ascendances, par unite de
75	!C surface de la maille (kg/m2/s)
76	!C p : pressions aux milieux des couches (Pa)
77	!C ph : pressions aux interfaces (Pa)
78	!C pi : (p/p_0)**kapa (adim)
79	!C dtime: increment temporel (s)
80	!c
81	!C Variables internes :
82	!c
83	!c rhow : masse volumique de la poche froide
84	!C rho : environment density at P levels
85	!C rhoh : environment density at Ph levels
86	!C te : environment temperature \| may change within
87	!C qe : environment humidity \| sub-time-stepping
88	!C the : environment potential temperature
89	!C thu : potential temperature in undisturbed area
90	!C tu : temperature in undisturbed area
91	!C qu : humidity in undisturbed area
92	!C dp_omgb: vertical gradient og LS omega
93	!C omgbw : wake average vertical omega
94	!C dp_omgbw: vertical gradient of omgbw
95	!C omgbdq : flux of Delta_q transported by LS omega
96	!C dth : potential temperature diff. wake-undist.
97	!C th1 : first pot. temp. for vertical advection (=thu)
98	!C th2 : second pot. temp. for vertical advection (=thw)
99	!C q1 : first humidity for vertical advection
100	!C q2 : second humidity for vertical advection
101	!C d_deltatw : terme de redistribution pour deltatw
102	!C d_deltaqw : terme de redistribution pour deltaqw
103	!C deltatw0 : deltatw initial
104	!C deltaqw0 : deltaqw initial
105	!C hw0 : hw initial
106	!C sigmaw0: sigmaw initial
107	!C amflux : horizontal mass flux through wake boundary
108	!C wdens_ref: initial number of wakes per unit area (3D) or per
109	!C unit length (2D), at the beginning of each time step
110	!C Tgw : 1 sur la période de onde de gravité
111	!c Cgw : vitesse de propagation de onde de gravité
112	!c LL : distance entre 2 poches
113
114	!c-------------------------------------------------------------------------
115	!c Déclaration de variables
116	!c-------------------------------------------------------------------------
117
118	#include "dimensions.h"
119	#include "YOMCST.h"
120	#include "cvthermo.h"
121	#include "iniprint.h"
122
123	!c Arguments en entree
124	!c--------------------
125
126	REAL, dimension(klon,klev) :: p, pi
127	REAL, dimension(klon,klev+1) :: ph, omgb
128	REAL dtime
129	REAL, dimension(klon,klev) :: te0,qe0
130	REAL, dimension(klon,klev) :: dtdwn, dqdwn
131	REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
132	REAL, dimension(klon,klev) :: udtPBL,udqPBL
133	REAL, dimension(klon,klev) :: amdwn, amup
134	REAL, dimension(klon,klev) :: dta, dqa
135	REAL, dimension(klon) :: sigd_con
136
137	!c Sorties
138	!c--------
139
140	REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
141	REAL, dimension(klon,klev) :: tu, qu
142	REAL, dimension(klon,klev) :: dtls, dqls
143	REAL, dimension(klon,klev) :: dtKE, dqKE
144	REAL, dimension(klon,klev) :: dtPBL, dqPBL
145	REAL, dimension(klon,klev) :: spread
146	REAL, dimension(klon,klev) :: d_deltatgw
147	REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
148	REAL, dimension(klon,klev+1) :: omgbdth, omg
149	REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
150	REAL, dimension(klon,klev) :: d_deltat_gw
151	REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
152	REAL, dimension(klon) :: wdens
153	INTEGER, dimension(klon) :: ktopw
154
155	!c Variables internes
156	!c-------------------
157
158	!c Variables à fixer
159	REAL ALON
160	REAL coefgw
161	REAL :: wdens_ref
162	REAL stark
163	REAL alpk
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	!c 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	!c 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	!c 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	!c 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	!c 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	REAL Crep_upper, Crep_sol
242
243	REAL, DIMENSION(klon,klev) :: ppi
244
245	!ccc nrlmd
246	real, dimension(klon) :: death_rate,nat_rate
247	real, dimension(klon,klev) :: entr
248	real, dimension(klon,klev) :: detr
249
250	!C-------------------------------------------------------------------------
251	!c Initialisations
252	!c-------------------------------------------------------------------------
253
254	!c print*, 'wake initialisations'
255
256	!c Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
257	!c-------------------------------------------------------------------------
258
259	DATA wapecut,sigmad, hwmin /5.,.02,10./
260	!ccc nrlmd
261	DATA sigmaw_max /0.4/
262	DATA dens_rate /0.1/
263	!ccc
264	!C Longueur de maille (en m)
265	!c-------------------------------------------------------------------------
266
267	!c ALON = 3.e5
268	ALON = 1.e6
269
270
271	!C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
272	!c
273	!c coefgw : Coefficient pour les ondes de gravité
274	!c stark : Coefficient k dans Cstar=ksqrt(2WAPE)
275	!c wdens : Densité de poche froide par maille
276	!c-------------------------------------------------------------------------
277
278	!ccc nrlmd coefgw=10
279	!c coefgw=1
280	!c wdens0 = 1.0/(alon**2)
281	!ccc nrlmd wdens = 1.0/(alon**2)
282	!ccc nrlmd stark = 0.50
283	!cCRtest
284	!ccc nrlmd alpk=0.1
285	!c alpk = 1.0
286	!c alpk = 0.5
287	!c alpk = 0.05
288	!c
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	!ccc nrlmd Lecture du fichier wake_param.data
297	OPEN(99,file='wake_param.data',status='old', &
298	& 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	!c
307	!c Initialisation de toutes des densites a wdens_ref.
308	!c Les densites peuvent evoluer si les poches debordent
309	!c (voir au tout debut de la boucle sur les substeps)
310	wdens = wdens_ref
311	!c
312	!c print*,'stark',stark
313	!c print*,'alpk',alpk
314	!c print*,'wdens',wdens
315	!c print*,'coefgw',coefgw
316	!ccc
317	!C Minimum value for \|T_wake - T_undist\|. Used for wake top definition
318	!c-------------------------------------------------------------------------
319
320	delta_t_min = 0.2
321
322	!C 1. - Save initial values and initialize tendencies
323	!C --------------------------------------------------
324
325	DO k=1,klev
326	DO i=1, klon
327	ppi(i,k)=pi(i,k)
328	deltatw0(i,k) = deltatw(i,k)
329	deltaqw0(i,k)= deltaqw(i,k)
330	te(i,k) = te0(i,k)
331	qe(i,k) = qe0(i,k)
332	dtls(i,k) = 0.
333	dqls(i,k) = 0.
334	d_deltat_gw(i,k)=0.
335	d_te(i,k) = 0.
336	d_qe(i,k) = 0.
337	d_deltatw(i,k) = 0.
338	d_deltaqw(i,k) = 0.
339	!IM 060508 beg
340	d_deltatw2(i,k)=0.
341	d_deltaqw2(i,k)=0.
342	!IM 060508 end
343	ENDDO
344	ENDDO
345	!c sigmaw1=sigmaw
346	!c IF (sigd_con.GT.sigmaw1) THEN
347	!c print*, 'sigmaw,sigd_con', sigmaw, sigd_con
348	!c ENDIF
349	DO i=1, klon
350	!cc sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
351	sigmaw(i) = amax1(sigmaw(i),sigmad)
352	sigmaw(i) = amin1(sigmaw(i),0.99)
353	sigmaw0(i) = sigmaw(i)
354	wape(i) = 0.
355	wape2(i) = 0.
356	d_sigmaw(i) = 0.
357	ktopw(i) = 0
358	ENDDO
359	!C
360	!C
361	!C 2. - Prognostic part
362	!C --------------------
363	!C
364	!C
365	!C 2.1 - Undisturbed area and Wake integrals
366	!C ---------------------------------------------------------
367
368	DO i=1, klon
369	z(i) = 0.
370	ktop(i)=0
371	kupper(i) = 0
372	sum_thu(i) = 0.
373	sum_tu(i) = 0.
374	sum_qu(i) = 0.
375	sum_thvu(i) = 0.
376	sum_dth(i) = 0.
377	sum_dq(i) = 0.
378	sum_rho(i) = 0.
379	sum_dtdwn(i) = 0.
380	sum_dqdwn(i) = 0.
381
382	av_thu(i) = 0.
383	av_tu(i) =0.
384	av_qu(i) =0.
385	av_thvu(i) = 0.
386	av_dth(i) = 0.
387	av_dq(i) = 0.
388	av_rho(i) =0.
389	av_dtdwn(i) =0.
390	av_dqdwn(i) = 0.
391	ENDDO
392	!c
393	!c Distance between wakes
394	DO i = 1,klon
395	LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
396	ENDDO
397	!C Potential temperatures and humidity
398	!c----------------------------------------------------------
399	DO k =1,klev
400	DO i=1, klon
401	! write(,)'wake 1',i,k,rd,te(i,k)
402	rho(i,k) = p(i,k)/(rd*te(i,k))
403	! write(,)'wake 2',rho(i,k)
404	IF(k .eq. 1) THEN
405	! write(,)'wake 3',i,k,rd,te(i,k)
406	rhoh(i,k) = ph(i,k)/(rd*te(i,k))
407	! write(,)'wake 4',i,k,rd,te(i,k)
408	zhh(i,k)=0
409	ELSE
410	! write(,)'wake 5',rd,(te(i,k)+te(i,k-1))
411	rhoh(i,k) = ph(i,k)2./(rd(te(i,k)+te(i,k-1)))
412	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
413	zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
414	ENDIF
415	! write(,)'wake 7',ppi(i,k)
416	the(i,k) = te(i,k)/ppi(i,k)
417	thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
418	tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
419	qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i)
420	! write(,)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
421	rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
422	dth(i,k) = deltatw(i,k)/ppi(i,k)
423	ENDDO
424	ENDDO
425
426	DO k = 1, klev-1
427	DO i=1, klon
428	IF(k.eq.1) THEN
429	N2(i,k)=0
430	ELSE
431	N2(i,k)=amax1(0.,-RG*2/the(i,k)rho(i,k)*(the(i,k+1)- &
432	& the(i,k-1))/(p(i,k+1)-p(i,k-1)))
433	ENDIF
434	ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2
435
436	Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)
437	Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)
438	ENDDO
439	ENDDO
440
441	DO i=1, klon
442	N2(i,klev)=0
443	ZH(i,klev)=0
444	Cgw(i,klev)=0
445	Tgw(i,klev)=0
446	ENDDO
447
448	!c Calcul de la masse volumique moyenne de la colonne (bdlmd)
449	!c-----------------------------------------------------------------
450
451	DO k=1,klev
452	DO i=1, klon
453	epaisseur1(i,k)=0.
454	epaisseur2(i,k)=0.
455	ENDDO
456	ENDDO
457
458	DO i=1, klon
459	epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
460	epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
461	rhow_moyen(i,1) = rhow(i,1)
462	ENDDO
463
464	DO k = 2, klev
465	DO i=1, klon
466	epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
467	epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
468	rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+ &
469	& rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
470	ENDDO
471	ENDDO
472
473	!C
474	!C Choose an integration bound well above wake top
475	!c-----------------------------------------------------------------
476	!c
477	!C Pupper = 50000. ! melting level
478	!c Pupper = 60000.
479	!c Pupper = 80000. ! essais pour case_e
480	DO i = 1,klon
481	Pupper(i) = 0.6*ph(i,1)
482	Pupper(i) = max(Pupper(i), 45000.)
483	!ccc Pupper(i) = 60000.
484	ENDDO
485
486	!C
487	!C Determine Wake top pressure (Ptop) from buoyancy integral
488	!C --------------------------------------------------------
489	!c
490	!c-1/ Pressure of the level where dth becomes less than delta_t_min.
491
492	DO i=1,klon
493	ptop_provis(i)=ph(i,1)
494	ENDDO
495	DO k= 2,klev
496	DO i=1,klon
497	!c
498	!IM v3JYG; ptop_provis(i).LT. ph(i,1)
499	!c
500	IF (dth(i,k) .GT. -delta_t_min .and. &
501	& dth(i,k-1).LT. -delta_t_min .and. &
502	& ptop_provis(i).EQ. ph(i,1)) THEN
503	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) &
504	& - (dth(i,k-1)+delta_t_min)*p(i,k)) / &
505	& (dth(i,k) - dth(i,k-1))
506	ENDIF
507	ENDDO
508	ENDDO
509
510	!c-2/ dth integral
511
512	DO i=1,klon
513	sum_dth(i) = 0.
514	dthmin(i) = -delta_t_min
515	z(i) = 0.
516	ENDDO
517
518	DO k = 1,klev
519	DO i=1,klon
520	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
521	IF (dz(i) .gt. 0) THEN
522	z(i) = z(i)+dz(i)
523	sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
524	dthmin(i) = amin1(dthmin(i),dth(i,k))
525	ENDIF
526	ENDDO
527	ENDDO
528
529	!c-3/ height of triangle with area= sum_dth and base = dthmin
530
531	DO i=1,klon
532	hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
533	hw0(i) = amax1(hwmin,hw0(i))
534	ENDDO
535
536	!c-4/ now, get Ptop
537
538	DO i=1,klon
539	z(i) = 0.
540	ptop(i) = ph(i,1)
541	ENDDO
542
543	DO k = 1,klev
544	DO i=1,klon
545	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
546	IF (dz(i) .gt. 0) THEN
547	z(i) = z(i)+dz(i)
548	ptop(i) = ph(i,k)-rho(i,k)rgdz(i)
549	ENDIF
550	ENDDO
551	ENDDO
552
553
554	!C-5/ Determination de ktop et kupper
555
556	DO k=klev,1,-1
557	DO i=1,klon
558	IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
559	IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
560	ENDDO
561	ENDDO
562
563	!c On evite kupper = 1 et kupper = klev
564	DO i=1,klon
565	kupper(i) = max(kupper(i),2)
566	kupper(i) = min(kupper(i),klev-1)
567	ENDDO
568
569
570	!c-6/ Correct ktop and ptop
571
572	DO i = 1,klon
573	ptop_new(i)=ptop(i)
574	ENDDO
575	DO k= klev,2,-1
576	DO i=1,klon
577	IF (k .LE. ktop(i) .and. &
578	& ptop_new(i) .EQ. ptop(i) .and. &
579	& dth(i,k) .GT. -delta_t_min .and. &
580	& dth(i,k-1).LT. -delta_t_min) THEN
581	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) &
582	& - (dth(i,k-1)+delta_t_min)*p(i,k)) / &
583	& (dth(i,k) - dth(i,k-1))
584	ENDIF
585	ENDDO
586	ENDDO
587
588	DO i=1,klon
589	ptop(i) = ptop_new(i)
590	ENDDO
591
592	DO k=klev,1,-1
593	DO i=1,klon
594	IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
595	ENDDO
596	ENDDO
597	!c
598	!c-5/ Set deltatw & deltaqw to 0 above kupper
599	!c
600	DO k = 1,klev
601	DO i=1,klon
602	IF (k.GE. kupper(i)) THEN
603	deltatw(i,k) = 0.
604	deltaqw(i,k) = 0.
605	ENDIF
606	ENDDO
607	ENDDO
608	!c
609	!C
610	!C Vertical gradient of LS omega
611	!C
612	DO k = 1,klev
613	DO i=1,klon
614	IF (k.LE. kupper(i)) THEN
615	dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
616	ENDIF
617	ENDDO
618	ENDDO
619	!C
620	!C Integrals (and wake top level number)
621	!C --------------------------------------
622	!C
623	!C Initialize sum_thvu to 1st level virt. pot. temp.
624
625	DO i=1,klon
626	z(i) = 1.
627	dz(i) = 1.
628	sum_thvu(i) = thu(i,1)(1.+epsqu(i,1))*dz(i)
629	sum_dth(i) = 0.
630	ENDDO
631
632	DO k = 1,klev
633	DO i=1,klon
634	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
635	IF (dz(i) .GT. 0) THEN
636	z(i) = z(i)+dz(i)
637	sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
638	sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
639	sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
640	sum_thvu(i) = sum_thvu(i) + thu(i,k)(1.+epsqu(i,k))*dz(i)
641	sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
642	sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
643	sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
644	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
645	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
646	ENDIF
647	ENDDO
648	ENDDO
649	!c
650	DO i=1,klon
651	hw0(i) = z(i)
652	ENDDO
653	!c
654	!C
655	!C 2.1 - WAPE and mean forcing computation
656	!C ---------------------------------------
657	!C
658	!C ---------------------------------------
659	!C
660	!C Means
661
662	DO i=1,klon
663	av_thu(i) = sum_thu(i)/hw0(i)
664	av_tu(i) = sum_tu(i)/hw0(i)
665	av_qu(i) = sum_qu(i)/hw0(i)
666	av_thvu(i) = sum_thvu(i)/hw0(i)
667	!c av_thve = sum_thve/hw0
668	av_dth(i) = sum_dth(i)/hw0(i)
669	av_dq(i) = sum_dq(i)/hw0(i)
670	av_rho(i) = sum_rho(i)/hw0(i)
671	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
672	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
673
674	wape(i) = - rghw0(i)(av_dth(i) &
675	& + eps(av_thu(i)av_dq(i)+av_dth(i)av_qu(i)+av_dth(i) &
676	& av_dq(i) ))/av_thvu(i)
677	ENDDO
678	!C
679	!C 2.2 Prognostic variable update
680	!C ------------------------------
681	!C
682	!C Filter out bad wakes
683
684	DO k = 1,klev
685	DO i=1,klon
686	IF ( wape(i) .LT. 0.) THEN
687	deltatw(i,k) = 0.
688	deltaqw(i,k) = 0.
689	dth(i,k) = 0.
690	ENDIF
691	ENDDO
692	ENDDO
693	!c
694	DO i=1,klon
695	IF ( wape(i) .LT. 0.) THEN
696	wape(i) = 0.
697	Cstar(i) = 0.
698	hw(i) = hwmin
699	sigmaw(i) = amax1(sigmad,sigd_con(i))
700	fip(i) = 0.
701	gwake(i) = .FALSE.
702	ELSE
703	Cstar(i) = starksqrt(2.wape(i))
704	gwake(i) = .TRUE.
705	ENDIF
706	ENDDO
707
708	!c
709	!c Check qx and qw positivity
710	!c --------------------------
711	DO i = 1,klon
712	q0_min(i)=min( (qe(i,1)-sigmaw(i)*deltaqw(i,1)), &
713	& (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)) )
714	ENDDO
715	DO k = 2,klev
716	DO i = 1,klon
717	q1_min(i)=min( (qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
718	& (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)) )
719	IF (q1_min(i).le.q0_min(i)) THEN
720	q0_min(i)=q1_min(i)
721	ENDIF
722	ENDDO
723	ENDDO
724	!c
725	DO i = 1,klon
726	OK_qx_qw(i) = q0_min(i) .GE. 0.
727	alpha(i) = 1.
728	ENDDO
729	!c
730	!CC -----------------------------------------------------------------
731	!C Sub-time-stepping
732	!C -----------------
733	!C
734	nsub=10
735	dtimesub=dtime/nsub
736	!c
737	!c------------------------------------------------------------
738	DO isubstep = 1,nsub
739	!c------------------------------------------------------------
740	!c
741	!c wk_adv is the logical flag enabling wake evolution in the time advance loop
742	DO i = 1,klon
743	wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1.
744	ENDDO
745	!c
746	!ccc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper --------
747	!ccc On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul ---
748	DO i=1,klon
749	!cc print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
750	!cc $ isubstep,wk_adv(i),cstar(i),wape(i)
751	IF (wk_adv(i) .AND. cstar(i).GT.0.01) THEN
752	omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i) &
753	& + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
754	wdens0 = ( sigmaw(i) / (4.3.14) ) &
755	& ( (1.-sigmaw(i)) * omg(i,kupper(i)+1) / &
756	& ( (ph(i,1)-pupper(i)) * cstar(i) ) ) **(2)
757	IF ( wdens(i) .LE. wdens0*1.1 ) THEN
758	wdens(i) = wdens0
759	ENDIF
760	!cc print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
761	!cc $ ,ph(i,1)-pupper(i)',
762	!cc $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
763	!cc $ ,ph(i,1)-pupper(i)
764	ENDIF
765	ENDDO
766
767	!ccc nrlmd
768
769	DO i=1,klon
770	IF (wk_adv(i)) THEN
771	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
772	sigmaw(i)=amin1(sigmaw(i),sigmaw_max)
773	ENDIF
774	ENDDO
775	DO i=1,klon
776	IF (wk_adv(i)) THEN
777	!ccc nrlmd Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4
778	!ccc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
779	IF (sigmaw(i).ge.sigmaw_max) THEN
780	death_rate(i)=gfl(i)*Cstar(i)/sigmaw(i)
781	ELSE
782	death_rate(i)=0.
783	END IF
784	d_sigmaw(i) = gfl(i)Cstar(i)dtimesub &
785	& - death_rate(i)sigmaw(i)dtimesub
786	!c $ - nat_rate(i)sigmaw(i)dtimesub
787	!cc print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
788	!cc $ death_rate(i),ktop(i),kupper(i)',
789	!cc $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
790	!cc $ death_rate(i),ktop(i),kupper(i)
791
792	!c sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
793	!c sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
794	!c wdens = wdens0/(10.*sigmaw)
795	!c sigmaw =max(sigmaw,sigd_con)
796	!c sigmaw =max(sigmaw,sigmad)
797	ENDIF
798	ENDDO
799	!C
800	!C
801	!c calcul de la difference de vitesse verticale poche - zone non perturbee
802	!IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
803	!IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
804	!IM 060208 au niveau k=1..?
805	DO k= 1,klev
806	DO i = 1,klon
807	if (wk_adv(i)) THEN !!! nrlmd
808	dp_deltomg(i,k)=0.
809	end if
810	ENDDO
811	ENDDO
812	DO k= 1,klev+1
813	DO i = 1,klon
814	if (wk_adv(i)) THEN !!! nrlmd
815	omg(i,k)=0.
816	end if
817	ENDDO
818	ENDDO
819	!c
820	DO i=1,klon
821	IF (wk_adv(i)) THEN
822	z(i)= 0.
823	omg(i,1) = 0.
824	dp_deltomg(i,1) = -(gfl(i)Cstar(i))/(sigmaw(i) (1-sigmaw(i)))
825	ENDIF
826	ENDDO
827	!c
828	DO k= 2,klev
829	DO i = 1,klon
830	IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
831	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
832	z(i) = z(i)+dz(i)
833	dp_deltomg(i,k)= dp_deltomg(i,1)
834	omg(i,k)= dp_deltomg(i,1)*z(i)
835	ENDIF
836	ENDDO
837	ENDDO
838	!c
839	DO i = 1,klon
840	IF (wk_adv(i)) THEN
841	dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
842	ztop(i) = z(i)+dztop(i)
843	omgtop(i)=dp_deltomg(i,1)*ztop(i)
844	ENDIF
845	ENDDO
846	!c
847	!c -----------------
848	!c From m/s to Pa/s
849	!c -----------------
850	!c
851	DO i=1,klon
852	IF (wk_adv(i)) THEN
853	omgtop(i) = -rho(i,ktop(i))rgomgtop(i)
854	dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
855	ENDIF
856	ENDDO
857	!c
858	DO k= 1,klev
859	DO i = 1,klon
860	IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
861	omg(i,k) = - rho(i,k)rgomg(i,k)
862	dp_deltomg(i,k) = dp_deltomg(i,1)
863	ENDIF
864	ENDDO
865	ENDDO
866	!c
867	!c raccordement lineaire de omg de ptop a pupper
868
869	DO i=1,klon
870	IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN
871	omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i) &
872	& + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
873	dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
874	& (ptop(i)-pupper(i))
875	ENDIF
876	ENDDO
877	!c
878	!cc DO i=1,klon
879	!cc print*,'Pente entre 0 et kupper (référence)'
880	!cc $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
881	!cc print*,'Pente entre ktop et kupper'
882	!cc $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
883	!cc ENDDO
884	!cc
885	DO k= 1,klev
886	DO i = 1,klon
887	IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
888	dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
889	omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
890	ENDIF
891	ENDDO
892	ENDDO
893	!ccc nrlmd
894	!cc DO i=1,klon
895	!cc print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
896	!cc END DO
897	!ccc
898	!c
899	!c
900	!c-- Compute wake average vertical velocity omgbw
901	!c
902	!c
903	DO k = 1,klev+1
904	DO i=1,klon
905	IF ( wk_adv(i)) THEN
906	omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
907	ENDIF
908	ENDDO
909	ENDDO
910	!c-- and its vertical gradient dp_omgbw
911	!c
912	DO k = 1,klev
913	DO i=1,klon
914	IF ( wk_adv(i)) THEN
915	dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
916	ENDIF
917	ENDDO
918	ENDDO
919	!C
920	!c-- Upstream coefficients for omgb velocity
921	!c-- (alpha_up(k) is the coefficient of the value at level k)
922	!c-- (1-alpha_up(k) is the coefficient of the value at level k-1)
923	DO k = 1,klev
924	DO i=1,klon
925	IF ( wk_adv(i)) THEN
926	alpha_up(i,k) = 0.
927	IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
928	ENDIF
929	ENDDO
930	ENDDO
931
932	!c Matrix expressing [The,deltatw] from [Th1,Th2]
933
934	DO i=1,klon
935	IF ( wk_adv(i)) THEN
936	RRe1(i) = 1.-sigmaw(i)
937	RRe2(i) = sigmaw(i)
938	ENDIF
939	ENDDO
940	RRd1 = -1.
941	RRd2 = 1.
942	!c
943	!c-- Get [Th1,Th2], dth and [q1,q2]
944	!c
945	DO k= 1,klev
946	DO i = 1,klon
947	IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
948	dth(i,k) = deltatw(i,k)/ppi(i,k)
949	Th1(i,k) = the(i,k) - sigmaw(i) *dth(i,k) ! undisturbed area
950	Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k) ! wake
951	q1(i,k) = qe(i,k) - sigmaw(i) *deltaqw(i,k) ! undisturbed area
952	q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
953	ENDIF
954	ENDDO
955	ENDDO
956
957	DO i=1,klon
958	if (wk_adv(i)) then !!! nrlmd
959	D_Th1(i,1) = 0.
960	D_Th2(i,1) = 0.
961	D_dth(i,1) = 0.
962	D_q1(i,1) = 0.
963	D_q2(i,1) = 0.
964	D_dq(i,1) = 0.
965	end if
966	ENDDO
967
968	DO k= 2,klev
969	DO i = 1,klon
970	IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
971	D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
972	D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
973	D_dth(i,k) = dth(i,k-1)-dth(i,k)
974	D_q1(i,k) = q1(i,k-1)-q1(i,k)
975	D_q2(i,k) = q2(i,k-1)-q2(i,k)
976	D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
977	ENDIF
978	ENDDO
979	ENDDO
980
981	DO i=1,klon
982	IF( wk_adv(i)) THEN
983	omgbdth(i,1) = 0.
984	omgbdq(i,1) = 0.
985	ENDIF
986	ENDDO
987
988	DO k= 2,klev
989	DO i = 1,klon
990	IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN ! loop on interfaces
991	omgbdth(i,k) = omgb(i,k)*( dth(i,k-1) - dth(i,k))
992	omgbdq(i,k) = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
993	ENDIF
994	ENDDO
995	ENDDO
996	!c
997	!c-----------------------------------------------------------------
998	DO k= 1,klev
999	DO i = 1,klon
1000	IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1001	!c-----------------------------------------------------------------
1002	!c
1003	!c Compute redistribution (advective) term
1004	!c
1005	d_deltatw(i,k) = &
1006	& dtimesub/(Ph(i,k)-Ph(i,k+1))*( &
1007	& RRd1omg(i,k )sigmaw(i) *D_Th1(i,k) &
1008	& -RRd2omg(i,k+1)(1.-sigmaw(i))*D_Th2(i,k+1) &
1009	& -(1.-alpha_up(i,k))omgbdth(i,k) - alpha_up(i,k+1) &
1010	& omgbdth(i,k+1))*ppi(i,k)
1011	!c print*,'d_deltatw=',d_deltatw(i,k)
1012	!c
1013	d_deltaqw(i,k) = &
1014	& dtimesub/(Ph(i,k)-Ph(i,k+1))*( &
1015	& RRd1omg(i,k )sigmaw(i) *D_q1(i,k) &
1016	& -RRd2omg(i,k+1)(1.-sigmaw(i))*D_q2(i,k+1) &
1017	& -(1.-alpha_up(i,k))omgbdq(i,k) - alpha_up(i,k+1) &
1018	& omgbdq(i,k+1))
1019	!c print*,'d_deltaqw=',d_deltaqw(i,k)
1020	!c
1021	!c and increment large scale tendencies
1022	!c
1023
1024	!c
1025	!C
1026	!CC -----------------------------------------------------------------
1027	d_te(i,k) = dtimesub*( &
1028	& ( RRe1(i)omg(i,k )sigmaw(i) *D_Th1(i,k) &
1029	& -RRe2(i)omg(i,k+1)(1.-sigmaw(i))*D_Th2(i,k+1) ) &
1030	& /(Ph(i,k)-Ph(i,k+1)) &
1031	!ccc nrlmd $ -sigmaw(i)(1.-sigmaw(i))dth(i,k)*dp_deltomg(i,k) &
1032	& -sigmaw(i)(1.-sigmaw(i))dth(i,k) &
1033	& *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) &
1034	!ccc &
1035	& )*ppi(i,k)
1036	!c
1037	d_qe(i,k) = dtimesub*( &
1038	& ( RRe1(i)omg(i,k )sigmaw(i) *D_q1(i,k) &
1039	& -RRe2(i)omg(i,k+1)(1.-sigmaw(i))*D_q2(i,k+1) ) &
1040	& /(Ph(i,k)-Ph(i,k+1)) &
1041	!ccc nrlmd $ -sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*dp_deltomg(i,k) &
1042	& -sigmaw(i)(1.-sigmaw(i))deltaqw(i,k) &
1043	& *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) &
1044	!ccc &
1045	& )
1046	!ccc nrlmd
1047	ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN
1048	d_te(i,k) = dtimesub*( &
1049	& ( RRe1(i)omg(i,k )sigmaw(i) *D_Th1(i,k) &
1050	& /(Ph(i,k)-Ph(i,k+1))) &
1051	& )*ppi(i,k)
1052
1053	d_qe(i,k) = dtimesub*( &
1054	& ( RRe1(i)omg(i,k )sigmaw(i) *D_q1(i,k) &
1055	& /(Ph(i,k)-Ph(i,k+1))) &
1056	& )
1057
1058	ENDIF
1059	!ccc
1060	ENDDO
1061	ENDDO
1062	!c------------------------------------------------------------------
1063	!C
1064	!C Increment state variables
1065
1066	DO k= 1,klev
1067	DO i = 1,klon
1068	!ccc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1069	IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1070	!ccc
1071
1072
1073	!c
1074	!c Coefficient de répartition
1075
1076	Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i)) &
1077	& -ph(i,1))
1078	Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)- &
1079	& ph(i,kupper(i)))
1080
1081
1082	!c Reintroduce compensating subsidence term.
1083
1084	!c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1085	!c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1086	!c . /(1-sigmaw)
1087	!c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1088	!c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1089	!c . /(1-sigmaw)
1090	!c
1091	!c dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
1092	!c dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
1093	!c . /(1-sigmaw)
1094	!c dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
1095	!c dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
1096	!c . /(1-sigmaw)
1097
1098	dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
1099	dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
1100	!c print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1101	!c
1102	dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))
1103	dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))
1104	!c print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
1105	!c
1106	!ccc nrlmd Prise en compte du taux de mortalité
1107	!ccc Définitions de entr, detr
1108	detr(i,k)=0.
1109
1110	entr(i,k)=detr(i,k)+gfl(i)*cstar(i)+ &
1111	& sigmaw(i)(1.-sigmaw(i))dp_deltomg(i,k)
1112
1113	spread(i,k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1114	!ccc spread(i,k) = (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
1115	!ccc $ sigmaw(i)
1116
1117
1118	!c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
1119
1120	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1121	! & Tgw(i,k),deltatw(i,k)
1122	d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)deltatw(i,k) &
1123	& dtimesub
1124	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1125	ff(i)=d_deltatw(i,k)/dtimesub
1126
1127	!c Sans GW
1128	!c
1129	!c deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
1130	!c
1131	!c GW formule 1
1132	!c
1133	!c deltatw(k) = deltatw(k)+dtimesub*
1134	!c $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
1135	!c
1136	!c GW formule 2
1137
1138	IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
1139	d_deltatw(i,k) = dtimesub* &
1140	& (ff(i)+dtKE(i,k)+dtPBL(i,k) &
1141	!ccc $ -spread(i,k)*deltatw(i,k) &
1142	& - entr(i,k)*deltatw(i,k)/sigmaw(i) &
1143	& - (death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k) &
1144	& / (1.-sigmaw(i)) &
1145	!ccc &
1146	& -Tgw(i,k)*deltatw(i,k))
1147	ELSE
1148	d_deltatw(i,k) = 1/Tgw(i,k)(1-exp(-dtimesub &
1149	& Tgw(i,k)))* &
1150	& (ff(i)+dtKE(i,k)+dtPBL(i,k) &
1151	!ccc $ -spread(i,k)*deltatw(i,k) &
1152	& - entr(i,k)*deltatw(i,k)/sigmaw(i) &
1153	& - (death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k) &
1154	& / (1.-sigmaw(i)) &
1155	!ccc &
1156	& -Tgw(i,k)*deltatw(i,k))
1157	ENDIF
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) &
1164	!ccc $ -spread(i,k)*deltaqw(i,k)) &
1165	& - entr(i,k)*deltaqw(i,k)/sigmaw(i) &
1166	& - (death_rate(i)sigmaw(i)+detr(i,k))deltaqw(i,k) &
1167	& /(1.-sigmaw(i)))
1168	!ccc
1169
1170	!ccc nrlmd
1171	!ccc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1172	!ccc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1173	!ccc
1174	ENDIF
1175	ENDDO
1176	ENDDO
1177
1178	!C
1179	!C Scale tendencies so that water vapour remains positive in w and x.
1180	!C
1181	call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw, &
1182	& d_deltaqw,sigmaw,d_sigmaw,alpha)
1183	!c
1184	!ccc nrlmd
1185	!cc print*,'alpha'
1186	!cc do i=1,klon
1187	!cc print*,alpha(i)
1188	!cc end do
1189	!ccc
1190	DO k = 1,klev
1191	DO i = 1,klon
1192	IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1193	d_te(i,k)=alpha(i)*d_te(i,k)
1194	d_qe(i,k)=alpha(i)*d_qe(i,k)
1195	d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)
1196	d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)
1197	d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)
1198	ENDIF
1199	ENDDO
1200	ENDDO
1201	DO i = 1,klon
1202	IF( wk_adv(i)) THEN
1203	d_sigmaw(i)=alpha(i)*d_sigmaw(i)
1204	ENDIF
1205	ENDDO
1206
1207	!C Update large scale variables and wake variables
1208	!IM 060208 manque DO i + remplace DO k=1,kupper(i)
1209	!IM 060208 DO k = 1,kupper(i)
1210	DO k= 1,klev
1211	DO i = 1,klon
1212	IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1213	dtls(i,k)=dtls(i,k)+d_te(i,k)
1214	dqls(i,k)=dqls(i,k)+d_qe(i,k)
1215	!ccc nrlmd
1216	d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1217	d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1218	!ccc
1219	ENDIF
1220	ENDDO
1221	ENDDO
1222	DO k= 1,klev
1223	DO i = 1,klon
1224	IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1225	te(i,k) = te0(i,k) + dtls(i,k)
1226	qe(i,k) = qe0(i,k) + dqls(i,k)
1227	the(i,k) = te(i,k)/ppi(i,k)
1228	deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)
1229	deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)
1230	dth(i,k) = deltatw(i,k)/ppi(i,k)
1231	!cc print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
1232	!cc $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1233	ENDIF
1234	ENDDO
1235	ENDDO
1236	DO i = 1,klon
1237	IF( wk_adv(i)) THEN
1238	sigmaw(i) = sigmaw(i)+d_sigmaw(i)
1239	ENDIF
1240	ENDDO
1241	!c
1242	!C
1243	!c Determine Ptop from buoyancy integral
1244	!c ---------------------------------------
1245	!c
1246	!c- 1/ Pressure of the level where dth changes sign.
1247	!c
1248	DO i=1,klon
1249	IF ( wk_adv(i)) THEN
1250	Ptop_provis(i)=ph(i,1)
1251	ENDIF
1252	ENDDO
1253	!c
1254	DO k= 2,klev
1255	DO i=1,klon
1256	IF ( wk_adv(i) .AND. &
1257	& Ptop_provis(i) .EQ. ph(i,1) .AND. &
1258	& dth(i,k) .GT. -delta_t_min .and. &
1259	& dth(i,k-1).LT. -delta_t_min) THEN
1260	Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) &
1261	& - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) &
1262	& - dth(i,k-1))
1263	ENDIF
1264	ENDDO
1265	ENDDO
1266	!c
1267	!c- 2/ dth integral
1268	!c
1269	DO i=1,klon
1270	if (wk_adv(i)) then !!! nrlmd
1271	sum_dth(i) = 0.
1272	dthmin(i) = -delta_t_min
1273	z(i) = 0.
1274	end if
1275	ENDDO
1276
1277	DO k = 1,klev
1278	DO i=1,klon
1279	IF ( wk_adv(i)) THEN
1280	dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
1281	IF (dz(i) .gt. 0) THEN
1282	z(i) = z(i)+dz(i)
1283	sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
1284	dthmin(i) = amin1(dthmin(i),dth(i,k))
1285	ENDIF
1286	ENDIF
1287	ENDDO
1288	ENDDO
1289	!c
1290	!c- 3/ height of triangle with area= sum_dth and base = dthmin
1291
1292	DO i=1,klon
1293	IF ( wk_adv(i)) THEN
1294	hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
1295	hw(i) = amax1(hwmin,hw(i))
1296	ENDIF
1297	ENDDO
1298	!c
1299	!c- 4/ now, get Ptop
1300	!c
1301	DO i=1,klon
1302	if (wk_adv(i)) then !!! nrlmd
1303	ktop(i) = 0
1304	z(i)=0.
1305	end if
1306	ENDDO
1307	!c
1308	DO k = 1,klev
1309	DO i=1,klon
1310	IF ( wk_adv(i)) THEN
1311	dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
1312	IF (dz(i) .gt. 0) THEN
1313	z(i) = z(i)+dz(i)
1314	Ptop(i) = Ph(i,k)-rho(i,k)rgdz(i)
1315	ktop(i) = k
1316	ENDIF
1317	ENDIF
1318	ENDDO
1319	ENDDO
1320	!c
1321	!c 4.5/Correct ktop and ptop
1322	!c
1323	DO i=1,klon
1324	IF ( wk_adv(i)) THEN
1325	Ptop_new(i)=ptop(i)
1326	ENDIF
1327	ENDDO
1328	!c
1329	DO k= klev,2,-1
1330	DO i=1,klon
1331	!IM v3JYG; IF (k .GE. ktop(i)
1332	IF ( wk_adv(i) .AND. &
1333	& k .LE. ktop(i) .AND. &
1334	& ptop_new(i) .EQ. ptop(i) .AND. &
1335	& dth(i,k) .GT. -delta_t_min .and. &
1336	& dth(i,k-1).LT. -delta_t_min) THEN
1337	Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) &
1338	& - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) &
1339	& - dth(i,k-1))
1340	ENDIF
1341	ENDDO
1342	ENDDO
1343	!c
1344	!c
1345	DO i=1,klon
1346	IF ( wk_adv(i)) THEN
1347	ptop(i) = ptop_new(i)
1348	ENDIF
1349	ENDDO
1350
1351	DO k=klev,1,-1
1352	DO i=1,klon
1353	if (wk_adv(i)) then !!! nrlmd
1354	IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
1355	end if
1356	ENDDO
1357	ENDDO
1358	!c
1359	!c 5/ Set deltatw & deltaqw to 0 above kupper
1360	!c
1361	DO k = 1,klev
1362	DO i=1,klon
1363	IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
1364	deltatw(i,k) = 0.
1365	deltaqw(i,k) = 0.
1366	ENDIF
1367	ENDDO
1368	ENDDO
1369	!c
1370	!C
1371	!c-------------Cstar computation---------------------------------
1372	DO i=1, klon
1373	if (wk_adv(i)) then !!! nrlmd
1374	sum_thu(i) = 0.
1375	sum_tu(i) = 0.
1376	sum_qu(i) = 0.
1377	sum_thvu(i) = 0.
1378	sum_dth(i) = 0.
1379	sum_dq(i) = 0.
1380	sum_rho(i) = 0.
1381	sum_dtdwn(i) = 0.
1382	sum_dqdwn(i) = 0.
1383
1384	av_thu(i) = 0.
1385	av_tu(i) =0.
1386	av_qu(i) =0.
1387	av_thvu(i) = 0.
1388	av_dth(i) = 0.
1389	av_dq(i) = 0.
1390	av_rho(i) =0.
1391	av_dtdwn(i) =0.
1392	av_dqdwn(i) = 0.
1393	end if
1394	ENDDO
1395	!C
1396	!C Integrals (and wake top level number)
1397	!C --------------------------------------
1398	!C
1399	!C Initialize sum_thvu to 1st level virt. pot. temp.
1400
1401	DO i=1,klon
1402	if (wk_adv(i)) then !!! nrlmd
1403	z(i) = 1.
1404	dz(i) = 1.
1405	sum_thvu(i) = thu(i,1)(1.+epsqu(i,1))*dz(i)
1406	sum_dth(i) = 0.
1407	end if
1408	ENDDO
1409
1410	DO k = 1,klev
1411	DO i=1,klon
1412	if (wk_adv(i)) then !!! nrlmd
1413	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1414	IF (dz(i) .GT. 0) THEN
1415	z(i) = z(i)+dz(i)
1416	sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
1417	sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
1418	sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
1419	sum_thvu(i) = sum_thvu(i) + thu(i,k)(1.+epsqu(i,k))*dz(i)
1420	sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
1421	sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
1422	sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
1423	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
1424	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
1425	ENDIF
1426	end if
1427	ENDDO
1428	ENDDO
1429	!c
1430	DO i=1,klon
1431	if (wk_adv(i)) then !!! nrlmd
1432	hw0(i) = z(i)
1433	end if
1434	ENDDO
1435	!c
1436	!C
1437	!C - WAPE and mean forcing computation
1438	!C ---------------------------------------
1439	!C
1440	!C ---------------------------------------
1441	!C
1442	!C Means
1443
1444	DO i=1,klon
1445	if (wk_adv(i)) then !!! nrlmd
1446	av_thu(i) = sum_thu(i)/hw0(i)
1447	av_tu(i) = sum_tu(i)/hw0(i)
1448	av_qu(i) = sum_qu(i)/hw0(i)
1449	av_thvu(i) = sum_thvu(i)/hw0(i)
1450	av_dth(i) = sum_dth(i)/hw0(i)
1451	av_dq(i) = sum_dq(i)/hw0(i)
1452	av_rho(i) = sum_rho(i)/hw0(i)
1453	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1454	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1455	!c
1456	wape(i) = - rghw0(i)(av_dth(i) &
1457	& + eps(av_thu(i)av_dq(i)+av_dth(i)av_qu(i)+av_dth(i) &
1458	& av_dq(i) ))/av_thvu(i)
1459	end if
1460	ENDDO
1461	!C
1462	!C Filter out bad wakes
1463
1464	DO k = 1,klev
1465	DO i=1,klon
1466	if (wk_adv(i)) then !!! nrlmd
1467	IF ( wape(i) .LT. 0.) THEN
1468	deltatw(i,k) = 0.
1469	deltaqw(i,k) = 0.
1470	dth(i,k) = 0.
1471	ENDIF
1472	end if
1473	ENDDO
1474	ENDDO
1475	!c
1476	DO i=1,klon
1477	if (wk_adv(i)) then !!! nrlmd
1478	IF ( wape(i) .LT. 0.) THEN
1479	wape(i) = 0.
1480	Cstar(i) = 0.
1481	hw(i) = hwmin
1482	sigmaw(i) = max(sigmad,sigd_con(i))
1483	fip(i) = 0.
1484	gwake(i) = .FALSE.
1485	ELSE
1486	Cstar(i) = starksqrt(2.wape(i))
1487	gwake(i) = .TRUE.
1488	ENDIF
1489	end if
1490	ENDDO
1491
1492	ENDDO ! end sub-timestep loop
1493	!C
1494	!C -----------------------------------------------------------------
1495	!c Get back to tendencies per second
1496	!c
1497	DO k = 1,klev
1498	DO i=1,klon
1499
1500	!ccc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1501	IF ( OK_qx_qw(i) .AND. k .LE. kupper(i)) THEN
1502	!ccc
1503	dtls(i,k) = dtls(i,k)/dtime
1504	dqls(i,k) = dqls(i,k)/dtime
1505	d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
1506	d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
1507	d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
1508	!cc print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
1509	!cc $ ,death_rate(i)*sigmaw(i)
1510	ENDIF
1511	ENDDO
1512	ENDDO
1513
1514	!c
1515	!c----------------------------------------------------------
1516	!c Determine wake final state; recompute wape, cstar, ktop;
1517	!c filter out bad wakes.
1518	!c----------------------------------------------------------
1519	!c
1520	!C 2.1 - Undisturbed area and Wake integrals
1521	!C ---------------------------------------------------------
1522
1523	DO i=1,klon
1524	!ccc nrlmd if (wk_adv(i)) then !!! nrlmd
1525	if (OK_qx_qw(i)) then
1526	!ccc
1527	z(i) = 0.
1528	sum_thu(i) = 0.
1529	sum_tu(i) = 0.
1530	sum_qu(i) = 0.
1531	sum_thvu(i) = 0.
1532	sum_dth(i) = 0.
1533	sum_dq(i) = 0.
1534	sum_rho(i) = 0.
1535	sum_dtdwn(i) = 0.
1536	sum_dqdwn(i) = 0.
1537
1538	av_thu(i) = 0.
1539	av_tu(i) =0.
1540	av_qu(i) =0.
1541	av_thvu(i) = 0.
1542	av_dth(i) = 0.
1543	av_dq(i) = 0.
1544	av_rho(i) =0.
1545	av_dtdwn(i) =0.
1546	av_dqdwn(i) = 0.
1547	end if
1548	ENDDO
1549	!C Potential temperatures and humidity
1550	!c----------------------------------------------------------
1551
1552	DO k =1,klev
1553	DO i=1,klon
1554	!ccc nrlmd IF ( wk_adv(i)) THEN
1555	if (OK_qx_qw(i)) then
1556	!ccc
1557	rho(i,k) = p(i,k)/(rd*te(i,k))
1558	IF(k .eq. 1) THEN
1559	rhoh(i,k) = ph(i,k)/(rd*te(i,k))
1560	zhh(i,k)=0
1561	ELSE
1562	rhoh(i,k) = ph(i,k)2./(rd(te(i,k)+te(i,k-1)))
1563	zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
1564	ENDIF
1565	the(i,k) = te(i,k)/ppi(i,k)
1566	thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
1567	tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
1568	qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i)
1569	rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
1570	dth(i,k) = deltatw(i,k)/ppi(i,k)
1571	ENDIF
1572	ENDDO
1573	ENDDO
1574
1575	!C Integrals (and wake top level number)
1576	!C -----------------------------------------------------------
1577
1578	!C Initialize sum_thvu to 1st level virt. pot. temp.
1579
1580	DO i=1,klon
1581	!ccc nrlmd IF ( wk_adv(i)) THEN
1582	if (OK_qx_qw(i)) then
1583	!ccc
1584	z(i) = 1.
1585	dz(i) = 1.
1586	sum_thvu(i) = thu(i,1)(1.+epsqu(i,1))*dz(i)
1587	sum_dth(i) = 0.
1588	ENDIF
1589	ENDDO
1590
1591	DO k = 1,klev
1592	DO i=1,klon
1593	!ccc nrlmd IF ( wk_adv(i)) THEN
1594	if (OK_qx_qw(i)) then
1595	!ccc
1596	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1597	IF (dz(i) .GT. 0) THEN
1598	z(i) = z(i)+dz(i)
1599	sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
1600	sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
1601	sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
1602	sum_thvu(i) = sum_thvu(i) + thu(i,k)(1.+epsqu(i,k))*dz(i)
1603	sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
1604	sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
1605	sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
1606	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
1607	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
1608	ENDIF
1609	ENDIF
1610	ENDDO
1611	ENDDO
1612	!c
1613	DO i=1,klon
1614	!ccc nrlmd IF ( wk_adv(i)) THEN
1615	if (OK_qx_qw(i)) then
1616	!ccc
1617	hw0(i) = z(i)
1618	ENDIF
1619	ENDDO
1620	!c
1621	!C - WAPE and mean forcing computation
1622	!C-------------------------------------------------------------
1623
1624	!C Means
1625
1626	DO i=1, klon
1627	!ccc nrlmd IF ( wk_adv(i)) THEN
1628	if (OK_qx_qw(i)) then
1629	!ccc
1630	av_thu(i) = sum_thu(i)/hw0(i)
1631	av_tu(i) = sum_tu(i)/hw0(i)
1632	av_qu(i) = sum_qu(i)/hw0(i)
1633	av_thvu(i) = sum_thvu(i)/hw0(i)
1634	av_dth(i) = sum_dth(i)/hw0(i)
1635	av_dq(i) = sum_dq(i)/hw0(i)
1636	av_rho(i) = sum_rho(i)/hw0(i)
1637	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1638	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1639
1640	wape2(i) = - rghw0(i)(av_dth(i) &
1641	& + eps(av_thu(i)av_dq(i)+av_dth(i)*av_qu(i) &
1642	& + av_dth(i)*av_dq(i) ))/av_thvu(i)
1643	ENDIF
1644	ENDDO
1645
1646	!C Prognostic variable update
1647	!C ------------------------------------------------------------
1648
1649	!C Filter out bad wakes
1650	!c
1651	DO k = 1,klev
1652	DO i=1,klon
1653	!ccc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
1654	if (OK_qx_qw(i) .AND. wape2(i) .LT. 0.) then
1655	!ccc
1656	deltatw(i,k) = 0.
1657	deltaqw(i,k) = 0.
1658	dth(i,k) = 0.
1659	ENDIF
1660	ENDDO
1661	ENDDO
1662	!c
1663
1664	DO i=1, klon
1665	!ccc nrlmd IF ( wk_adv(i)) THEN
1666	if (OK_qx_qw(i)) then
1667	!ccc
1668	IF ( wape2(i) .LT. 0.) THEN
1669	wape2(i) = 0.
1670	Cstar2(i) = 0.
1671	hw(i) = hwmin
1672	sigmaw(i) = amax1(sigmad,sigd_con(i))
1673	fip(i) = 0.
1674	gwake(i) = .FALSE.
1675	ELSE
1676	if(prt_level.ge.10) print*,'wape2>0'
1677	Cstar2(i) = starksqrt(2.wape2(i))
1678	gwake(i) = .TRUE.
1679	ENDIF
1680	ENDIF
1681	ENDDO
1682	!c
1683	DO i=1, klon
1684	!ccc nrlmd IF ( wk_adv(i)) THEN
1685	if (OK_qx_qw(i)) then
1686	!ccc
1687	ktopw(i) = ktop(i)
1688	ENDIF
1689	ENDDO
1690	!c
1691	DO i=1, klon
1692	!ccc nrlmd IF ( wk_adv(i)) THEN
1693	if (OK_qx_qw(i)) then
1694	!ccc
1695	IF (ktopw(i) .gt. 0 .and. gwake(i)) then
1696
1697	!Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
1698	!ccc heff = 600.
1699	!C Utilisation de la hauteur hw
1700	!cc heff = 0.7*hw
1701	heff(i) = hw(i)
1702
1703	FIP(i) = 0.5rho(i,ktopw(i))Cstar2(i)*3heff(i)2 &
1704	& sqrt(sigmaw(i)wdens(i)3.14)
1705	FIP(i) = alpk * FIP(i)
1706	!Cjyg2
1707	ELSE
1708	FIP(i) = 0.
1709	ENDIF
1710	ENDIF
1711	ENDDO
1712	!c
1713	!C Limitation de sigmaw
1714
1715	!ccc nrlmd
1716	!c DO i=1,klon
1717	!c IF (OK_qx_qw(i)) THEN
1718	!c IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
1719	!c ENDIF
1720	!c ENDDO
1721	!ccc
1722	DO k = 1,klev
1723	DO i=1, klon
1724
1725	!ccc nrlmd On maintient désormais constant sigmaw en régime permanent
1726	!ccc IF ((sigmaw(i).GT.sigmaw_max).or.
1727	IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. &
1728	& (ktopw(i).le.2) .OR. &
1729	& .not. OK_qx_qw(i) ) THEN
1730	!ccc
1731	dtls(i,k) = 0.
1732	dqls(i,k) = 0.
1733	deltatw(i,k) = 0.
1734	deltaqw(i,k) = 0.
1735	ENDIF
1736	ENDDO
1737	ENDDO
1738	!c
1739	!ccc nrlmd On maintient désormais constant sigmaw en régime permanent
1740	DO i=1, klon
1741	IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. &
1742	& (ktopw(i).le.2) .OR. &
1743	& .not. OK_qx_qw(i) ) THEN
1744	wape(i) = 0.
1745	cstar(i)=0.
1746	hw(i) = hwmin
1747	sigmaw(i) = sigmad
1748	fip(i) = 0.
1749	ELSE
1750	wape(i) = wape2(i)
1751	cstar(i)=cstar2(i)
1752	ENDIF
1753	!cc print*,'wape wape2 ktopw OK_qx_qw =',
1754	!cc $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
1755	ENDDO
1756	!c
1757	!c
1758	RETURN
1759	END SUBROUTINE WAKE
1760
1761	SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe, &
1762	& deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
1763	!c------------------------------------------------------
1764	!cDtermination du coefficient alpha tel que les tendances
1765	!c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
1766	!c a une humidite positive dans la zone (x) et dans la zone (w).
1767	!c------------------------------------------------------
1768	!c
1769
1770	!c Input
1771	REAL qe(nlon,nl),d_qe(nlon,nl)
1772	REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)
1773	REAL sigmaw(nlon),d_sigmaw(nlon)
1774	LOGICAL wk_adv(nlon)
1775	INTEGER nl,nlon
1776	!c Output
1777	REAL alpha(nlon)
1778	!c Internal variables
1779	REAL zeta(nlon,nl)
1780	REAL alpha1(nlon)
1781	REAL x,a,b,c,discrim
1782	REAL epsilon
1783	! DATA epsilon/1.e-15/
1784	!c
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)).ge.0.) then
1789	zeta(i,k)=0.
1790	ELSE
1791	zeta(i,k)=1.
1792	END IF
1793	ENDIF
1794	ENDDO
1795	DO i = 1,nlon
1796	IF (wk_adv(i)) THEN
1797	x = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k) &
1798	& + d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k) &
1799	& - d_sigmaw(i)*(deltaqw(i,k)+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	!c print*, 'x, a, b, c, discrim', x, a, b, c, discrim
1806	IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap
1807	alpha1(i)=1.
1808	ELSE
1809	IF (x .GE. 0.) THEN
1810	alpha1(i)=1.
1811	ELSE
1812	IF (a .GT. 0.) THEN
1813	alpha1(i)=0.9min( (2.c)/(-b+sqrt(discrim)), &
1814	& (-b+sqrt(discrim))/(2.*a) )
1815	ELSE IF (a .eq. 0.) then
1816	alpha1(i)=0.9*(-c/b)
1817	ELSE
1818	!c print*,'a,b,c discrim',a,b,c discrim
1819	alpha1(i)=0.9max( (2.c)/(-b+sqrt(discrim)), &
1820	& (-b+sqrt(discrim))/(2.*a) )
1821	ENDIF
1822	ENDIF
1823	ENDIF
1824	alpha(i) = min(alpha(i),alpha1(i))
1825	ENDIF
1826	ENDDO
1827	ENDDO
1828	!
1829	return
1830	end SUBROUTINE wake_vec_modulation
1831
1832	SUBROUTINE WAKE_scal (p,ph,ppi,dtime,sigd_con &
1833	& ,te0,qe0,omgb &
1834	& ,dtdwn,dqdwn,amdwn,amup,dta,dqa &
1835	& ,wdtPBL,wdqPBL,udtPBL,udqPBL &
1836	& ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl &
1837	& ,dtls,dqls &
1838	& ,ktopw,omgbdth,dp_omgb,wdens &
1839	& ,tu,qu &
1840	& ,dtKE,dqKE &
1841	& ,dtPBL,dqPBL &
1842	& ,omg,dp_deltomg,spread &
1843	& ,Cstar,d_deltat_gw &
1844	& ,d_deltatw2,d_deltaqw2)
1845
1846	!***************************************************************
1847	!* *
1848	!* WAKE *
1849	!* retour a un Pupper fixe *
1850	!* *
1851	!* written by : GRANDPEIX Jean-Yves 09/03/2000 *
1852	!* modified by : ROEHRIG Romain 01/29/2007 *
1853	!***************************************************************
1854	!c
1855	USE dimphy
1856	IMPLICIT none
1857	!c============================================================================
1858	!C
1859	!C
1860	!C But : Decrire le comportement des poches froides apparaissant dans les
1861	!C grands systemes convectifs, et fournir l'energie disponible pour
1862	!C le declenchement de nouvelles colonnes convectives.
1863	!C
1864	!C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area
1865	!C deltaqw : ecart d'humidite wake-undisturbed area
1866	!C sigmaw : fraction d'aire occupee par la poche.
1867	!C
1868	!C Variable de sortie :
1869	!c
1870	!c wape : WAke Potential Energy
1871	!c fip : Front Incident Power (W/m2) - ALP
1872	!c gfl : Gust Front Length per unit area (m-1)
1873	!C dtls : large scale temperature tendency due to wake
1874	!C dqls : large scale humidity tendency due to wake
1875	!C hw : hauteur de la poche
1876	!C dp_omgb : vertical gradient of large scale omega
1877	!C omgbdth: flux of Delta_Theta transported by LS omega
1878	!C dtKE : differential heating (wake - unpertubed)
1879	!C dqKE : differential moistening (wake - unpertubed)
1880	!C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
1881	!C dp_deltomg : vertical gradient of omg (s-1)
1882	!C spread : spreading term in dt_wake and dq_wake
1883	!C deltatw : updated temperature difference (T_w-T_u).
1884	!C deltaqw : updated humidity difference (q_w-q_u).
1885	!C sigmaw : updated wake fractional area.
1886	!C d_deltat_gw : delta T tendency due to GW
1887	!c
1888	!C Variables d'entree :
1889	!c
1890	!c aire : aire de la maille
1891	!c te0 : temperature dans l'environnement (K)
1892	!C qe0 : humidite dans l'environnement (kg/kg)
1893	!C omgb : vitesse verticale moyenne sur la maille (Pa/s)
1894	!C dtdwn: source de chaleur due aux descentes (K/s)
1895	!C dqdwn: source d'humidite due aux descentes (kg/kg/s)
1896	!C dta : source de chaleur due courants satures et detrain (K/s)
1897	!C dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
1898	!C amdwn: flux de masse total des descentes, par unite de
1899	!C surface de la maille (kg/m2/s)
1900	!C amup : flux de masse total des ascendances, par unite de
1901	!C surface de la maille (kg/m2/s)
1902	!C p : pressions aux milieux des couches (Pa)
1903	!C ph : pressions aux interfaces (Pa)
1904	!C ppi : (p/p_0)**kapa (adim)
1905	!C dtime: increment temporel (s)
1906	!c
1907	!C Variables internes :
1908	!c
1909	!c rhow : masse volumique de la poche froide
1910	!C rho : environment density at P levels
1911	!C rhoh : environment density at Ph levels
1912	!C te : environment temperature \| may change within
1913	!C qe : environment humidity \| sub-time-stepping
1914	!C the : environment potential temperature
1915	!C thu : potential temperature in undisturbed area
1916	!C tu : temperature in undisturbed area
1917	!C qu : humidity in undisturbed area
1918	!C dp_omgb: vertical gradient og LS omega
1919	!C omgbw : wake average vertical omega
1920	!C dp_omgbw: vertical gradient of omgbw
1921	!C omgbdq : flux of Delta_q transported by LS omega
1922	!C dth : potential temperature diff. wake-undist.
1923	!C th1 : first pot. temp. for vertical advection (=thu)
1924	!C th2 : second pot. temp. for vertical advection (=thw)
1925	!C q1 : first humidity for vertical advection
1926	!C q2 : second humidity for vertical advection
1927	!C d_deltatw : terme de redistribution pour deltatw
1928	!C d_deltaqw : terme de redistribution pour deltaqw
1929	!C deltatw0 : deltatw initial
1930	!C deltaqw0 : deltaqw initial
1931	!C hw0 : hw initial
1932	!C sigmaw0: sigmaw initial
1933	!C amflux : horizontal mass flux through wake boundary
1934	!C wdens : number of wakes per unit area (3D) or per
1935	!C unit length (2D)
1936	!C Tgw : 1 sur la période de onde de gravité
1937	!c Cgw : vitesse de propagation de onde de gravité
1938	!c LL : distance entre 2 poches
1939
1940	!c-------------------------------------------------------------------------
1941	!c Déclaration de variables
1942	!c-------------------------------------------------------------------------
1943
1944	#include "dimensions.h"
1945	!cccc#include "dimphy.h"
1946	#include "YOMCST.h"
1947	#include "cvthermo.h"
1948	#include "iniprint.h"
1949
1950	!c Arguments en entree
1951	!c--------------------
1952
1953	REAL p(klev),ph(klev+1),ppi(klev)
1954	REAL dtime
1955	REAL te0(klev),qe0(klev)
1956	REAL omgb(klev+1)
1957	REAL dtdwn(klev), dqdwn(klev)
1958	REAL wdtPBL(klev),wdqPBL(klev)
1959	REAL udtPBL(klev),udqPBL(klev)
1960	REAL amdwn(klev), amup(klev)
1961	REAL dta(klev), dqa(klev)
1962	REAL sigd_con
1963
1964	!c Sorties
1965	!c--------
1966
1967	REAL deltatw(klev), deltaqw(klev), dth(klev)
1968	REAL tu(klev), qu(klev)
1969	REAL dtls(klev), dqls(klev)
1970	REAL dtKE(klev), dqKE(klev)
1971	REAL dtPBL(klev), dqPBL(klev)
1972	REAL spread(klev)
1973	REAL d_deltatgw(klev)
1974	REAL d_deltatw2(klev), d_deltaqw2(klev)
1975	REAL omgbdth(klev+1), omg(klev+1)
1976	REAL dp_omgb(klev), dp_deltomg(klev)
1977	REAL d_deltat_gw(klev)
1978	REAL hw, sigmaw, wape, fip, gfl, Cstar
1979	INTEGER ktopw
1980
1981	!c Variables internes
1982	!c-------------------
1983
1984	!c Variables à fixer
1985	REAL ALON
1986	REAL coefgw
1987	REAL wdens0, wdens
1988	REAL stark
1989	REAL alpk
1990	REAL delta_t_min
1991	REAL Pupper
1992	INTEGER nsub
1993	REAL dtimesub
1994	REAL sigmad, hwmin
1995
1996	!c Variables de sauvegarde
1997	REAL deltatw0(klev)
1998	REAL deltaqw0(klev)
1999	REAL te(klev), qe(klev)
2000	REAL sigmaw0, sigmaw1
2001
2002	!c Variables pour les GW
2003	REAL LL
2004	REAL N2(klev)
2005	REAL Cgw(klev)
2006	REAL Tgw(klev)
2007
2008	!c Variables liées au calcul de hw
2009	REAL ptop_provis, ptop, ptop_new
2010	REAL sum_dth
2011	REAL dthmin
2012	REAL z, dz, hw0
2013	INTEGER ktop, kupper
2014
2015	!c Autres variables internes
2016	INTEGER isubstep, k
2017
2018	REAL sum_thu, sum_tu, sum_qu,sum_thvu
2019	REAL sum_dq, sum_rho
2020	REAL sum_dtdwn, sum_dqdwn
2021	REAL av_thu, av_tu, av_qu, av_thvu
2022	REAL av_dth, av_dq, av_rho
2023	REAL av_dtdwn, av_dqdwn
2024
2025	REAL rho(klev), rhoh(klev+1), rhow(klev)
2026	REAL rhow_moyen(klev)
2027	REAL zh(klev), zhh(klev+1)
2028	REAL epaisseur1(klev), epaisseur2(klev)
2029
2030	REAL the(klev), thu(klev)
2031
2032	REAL d_deltatw(klev), d_deltaqw(klev)
2033
2034	REAL omgbw(klev+1), omgtop
2035	REAL dp_omgbw(klev)
2036	REAL ztop, dztop
2037	REAL alpha_up(klev)
2038
2039	REAL RRe1, RRe2, RRd1, RRd2
2040	REAL Th1(klev), Th2(klev), q1(klev), q2(klev)
2041	REAL D_Th1(klev), D_Th2(klev), D_dth(klev)
2042	REAL D_q1(klev), D_q2(klev), D_dq(klev)
2043	REAL omgbdq(klev)
2044
2045	REAL ff, gg
2046	REAL wape2, Cstar2, heff
2047
2048	REAL Crep(klev)
2049	REAL Crep_upper, Crep_sol
2050
2051	!C-------------------------------------------------------------------------
2052	!c Initialisations
2053	!c-------------------------------------------------------------------------
2054
2055	!c print*, 'wake initialisations'
2056
2057	!c Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
2058	!c-------------------------------------------------------------------------
2059
2060	DATA sigmad, hwmin /.02,10./
2061
2062	!C Longueur de maille (en m)
2063	!c-------------------------------------------------------------------------
2064
2065	!c ALON = 3.e5
2066	ALON = 1.e6
2067
2068
2069	!C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
2070	!c
2071	!c coefgw : Coefficient pour les ondes de gravité
2072	!c stark : Coefficient k dans Cstar=ksqrt(2WAPE)
2073	!c wdens : Densité de poche froide par maille
2074	!c-------------------------------------------------------------------------
2075
2076	coefgw=10
2077	!c coefgw=1
2078	!c wdens0 = 1.0/(alon**2)
2079	wdens = 1.0/(alon**2)
2080	stark = 0.50
2081	!cCRtest
2082	alpk=0.1
2083	!c alpk = 1.0
2084	!c alpk = 0.5
2085	!c alpk = 0.05
2086	Crep_upper=0.9
2087	Crep_sol=1.0
2088
2089
2090	!C Minimum value for \|T_wake - T_undist\|. Used for wake top definition
2091	!c-------------------------------------------------------------------------
2092
2093	delta_t_min = 0.2
2094
2095
2096	!C 1. - Save initial values and initialize tendencies
2097	!C --------------------------------------------------
2098
2099	DO k=1,klev
2100	deltatw0(k) = deltatw(k)
2101	deltaqw0(k)= deltaqw(k)
2102	te(k) = te0(k)
2103	qe(k) = qe0(k)
2104	dtls(k) = 0.
2105	dqls(k) = 0.
2106	d_deltat_gw(k)=0.
2107	d_deltatw2(k)=0.
2108	d_deltaqw2(k)=0.
2109	ENDDO
2110	!c sigmaw1=sigmaw
2111	!c IF (sigd_con.GT.sigmaw1) THEN
2112	!c print*, 'sigmaw,sigd_con', sigmaw, sigd_con
2113	!c ENDIF
2114	sigmaw = max(sigmaw,sigd_con)
2115	sigmaw = max(sigmaw,sigmad)
2116	sigmaw = min(sigmaw,0.99)
2117	sigmaw0 = sigmaw
2118	!c wdens=wdens0/(10.*sigmaw)
2119	!c IF (sigd_con.GT.sigmaw1) THEN
2120	!c print*, 'sigmaw1,sigd1', sigmaw, sigd_con
2121	!c ENDIF
2122
2123	!C 2. - Prognostic part
2124	!C =========================================================
2125
2126	!c print *, 'prognostic wake computation'
2127
2128
2129	!C 2.1 - Undisturbed area and Wake integrals
2130	!C ---------------------------------------------------------
2131
2132	z = 0.
2133	ktop=0
2134	kupper = 0
2135	sum_thu = 0.
2136	sum_tu = 0.
2137	sum_qu = 0.
2138	sum_thvu = 0.
2139	sum_dth = 0.
2140	sum_dq = 0.
2141	sum_rho = 0.
2142	sum_dtdwn = 0.
2143	sum_dqdwn = 0.
2144
2145	av_thu = 0.
2146	av_tu =0.
2147	av_qu =0.
2148	av_thvu = 0.
2149	av_dth = 0.
2150	av_dq = 0.
2151	av_rho =0.
2152	av_dtdwn =0.
2153	av_dqdwn = 0.
2154
2155	!C Potential temperatures and humidity
2156	!c----------------------------------------------------------
2157
2158	DO k =1,klev
2159	rho(k) = p(k)/(rd*te(k))
2160	IF(k .eq. 1) THEN
2161	rhoh(k) = ph(k)/(rd*te(k))
2162	zhh(k)=0
2163	ELSE
2164	rhoh(k) = ph(k)2./(rd(te(k)+te(k-1)))
2165	zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
2166	ENDIF
2167	the(k) = te(k)/ppi(k)
2168	thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
2169	tu(k) = te(k) - deltatw(k)*sigmaw
2170	qu(k) = qe(k) - deltaqw(k)*sigmaw
2171	rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
2172	dth(k) = deltatw(k)/ppi(k)
2173	LL = (1-sqrt(sigmaw))/sqrt(wdens)
2174	ENDDO
2175
2176	DO k = 1, klev-1
2177	IF(k.eq.1) THEN
2178	N2(k)=0
2179	ELSE
2180	N2(k)=max(0.,-RG*2/the(k)rho(k)*(the(k+1)-the(k-1)) &
2181	& /(p(k+1)-p(k-1)))
2182	ENDIF
2183	ZH(k)=(zhh(k)+zhh(k+1))/2
2184
2185	Cgw(k)=sqrt(N2(k))*ZH(k)
2186	Tgw(k)=coefgw*Cgw(k)/LL
2187	ENDDO
2188
2189	N2(klev)=0
2190	ZH(klev)=0
2191	Cgw(klev)=0
2192	Tgw(klev)=0
2193
2194	!c Calcul de la masse volumique moyenne de la colonne
2195	!c-----------------------------------------------------------------
2196
2197	DO k=1,klev
2198	epaisseur1(k)=0.
2199	epaisseur2(k)=0.
2200	ENDDO
2201
2202	epaisseur1(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
2203	epaisseur2(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
2204	rhow_moyen(1) = rhow(1)
2205
2206	DO k = 2, klev
2207	epaisseur1(k)= -(Ph(k+1)-Ph(k))/(rho(k)*rg) +1.
2208	epaisseur2(k)=epaisseur2(k-1)+epaisseur1(k)
2209	rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+ &
2210	& rhow(k)*epaisseur1(k))/epaisseur2(k)
2211	ENDDO
2212
2213
2214	!C Choose an integration bound well above wake top
2215	!c-----------------------------------------------------------------
2216
2217	!c Pupper = 50000. ! melting level
2218	Pupper = 60000.
2219	!c Pupper = 70000.
2220
2221
2222	!C Determine Wake top pressure (Ptop) from buoyancy integral
2223	!C-----------------------------------------------------------------
2224
2225	!c-1/ Pressure of the level where dth becomes less than delta_t_min.
2226
2227	Ptop_provis=ph(1)
2228	DO k= 2,klev
2229	IF (dth(k) .GT. -delta_t_min .and. &
2230	& dth(k-1).LT. -delta_t_min) THEN
2231	Ptop_provis = ((dth(k)+delta_t_min)*p(k-1) &
2232	& - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2233	GO TO 25
2234	ENDIF
2235	ENDDO
2236	25 CONTINUE
2237
2238	!c-2/ dth integral
2239
2240	sum_dth = 0.
2241	dthmin = -delta_t_min
2242	z = 0.
2243
2244	DO k = 1,klev
2245	dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
2246	IF (dz .le. 0) GO TO 40
2247	z = z+dz
2248	sum_dth = sum_dth + dth(k)*dz
2249	dthmin = min(dthmin,dth(k))
2250	ENDDO
2251	40 CONTINUE
2252
2253	!c-3/ height of triangle with area= sum_dth and base = dthmin
2254
2255	hw0 = 2.*sum_dth/min(dthmin,-0.5)
2256	hw0 = max(hwmin,hw0)
2257
2258	!c-4/ now, get Ptop
2259
2260	z = 0.
2261	ptop = ph(1)
2262
2263	DO k = 1,klev
2264	dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw0-z)
2265	IF (dz .le. 0) GO TO 45
2266	z = z+dz
2267	Ptop = Ph(k)-rho(k)rgdz
2268	ENDDO
2269	45 CONTINUE
2270
2271
2272	!C-5/ Determination de ktop et kupper
2273
2274	DO k=klev,1,-1
2275	IF (ph(k+1) .lt. ptop) ktop=k
2276	IF (ph(k+1) .lt. pupper) kupper=k
2277	ENDDO
2278
2279	!c-6/ Correct ktop and ptop
2280
2281	Ptop_new=ptop
2282	DO k= ktop,2,-1
2283	IF (dth(k) .GT. -delta_t_min .and. &
2284	& dth(k-1).LT. -delta_t_min) THEN
2285	Ptop_new = ((dth(k)+delta_t_min)*p(k-1) &
2286	& - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2287	GO TO 225
2288	ENDIF
2289	ENDDO
2290	225 CONTINUE
2291
2292	ptop = ptop_new
2293
2294	DO k=klev,1,-1
2295	IF (ph(k+1) .lt. ptop) ktop=k
2296	ENDDO
2297
2298	!c Set deltatw & deltaqw to 0 above kupper
2299	!c-----------------------------------------------------------
2300
2301	DO k = kupper,klev
2302	deltatw(k) = 0.
2303	deltaqw(k) = 0.
2304	ENDDO
2305
2306
2307	!C Vertical gradient of LS omega
2308	!C------------------------------------------------------------
2309
2310	DO k = 1,kupper
2311	dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k))
2312	ENDDO
2313
2314
2315	!C Integrals (and wake top level number)
2316	!C -----------------------------------------------------------
2317
2318	!C Initialize sum_thvu to 1st level virt. pot. temp.
2319
2320	z = 1.
2321	dz = 1.
2322	sum_thvu = thu(1)(1.+epsqu(1))*dz
2323	sum_dth = 0.
2324
2325	DO k = 1,klev
2326	dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
2327	IF (dz .LE. 0) GO TO 50
2328	z = z+dz
2329	sum_thu = sum_thu + thu(k)*dz
2330	sum_tu = sum_tu + tu(k)*dz
2331	sum_qu = sum_qu + qu(k)*dz
2332	sum_thvu = sum_thvu + thu(k)(1.+epsqu(k))*dz
2333	sum_dth = sum_dth + dth(k)*dz
2334	sum_dq = sum_dq + deltaqw(k)*dz
2335	sum_rho = sum_rho + rhow(k)*dz
2336	sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
2337	sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
2338	ENDDO
2339	50 CONTINUE
2340
2341	hw0 = z
2342
2343	!C 2.1 - WAPE and mean forcing computation
2344	!C-------------------------------------------------------------
2345
2346	!C Means
2347
2348	av_thu = sum_thu/hw0
2349	av_tu = sum_tu/hw0
2350	av_qu = sum_qu/hw0
2351	av_thvu = sum_thvu/hw0
2352	!c av_thve = sum_thve/hw0
2353	av_dth = sum_dth/hw0
2354	av_dq = sum_dq/hw0
2355	av_rho = sum_rho/hw0
2356	av_dtdwn = sum_dtdwn/hw0
2357	av_dqdwn = sum_dqdwn/hw0
2358
2359	wape = - rghw0(av_dth &
2360	& + eps(av_thuav_dq+av_dthav_qu+av_dthav_dq ))/av_thvu
2361
2362	!C 2.2 Prognostic variable update
2363	!C ------------------------------------------------------------
2364
2365	!C Filter out bad wakes
2366
2367	IF ( wape .LT. 0.) THEN
2368	if(prt_level.ge.10) print*,'wape<0'
2369	wape = 0.
2370	hw = hwmin
2371	sigmaw = max(sigmad,sigd_con)
2372	fip = 0.
2373	DO k = 1,klev
2374	deltatw(k) = 0.
2375	deltaqw(k) = 0.
2376	dth(k) = 0.
2377	ENDDO
2378	ELSE
2379	if(prt_level.ge.10) print*,'wape>0'
2380	Cstar = starksqrt(2.wape)
2381	ENDIF
2382
2383	!C------------------------------------------------------------------
2384	!C Sub-time-stepping
2385	!C------------------------------------------------------------------
2386
2387	!c nsub=36
2388	nsub=10
2389	dtimesub=dtime/nsub
2390
2391	!c------------------------------------------------------------
2392	DO isubstep = 1,nsub
2393	!c------------------------------------------------------------
2394
2395	!c print*,'---------------','substep=',isubstep,'-------------'
2396
2397	!c Evolution of sigmaw
2398
2399
2400	gfl = 2.sqrt(3.14wdens*sigmaw)
2401
2402	sigmaw =sigmaw + gflCstardtimesub
2403	sigmaw =min(sigmaw,0.99) !!!!!!!!
2404	!c wdens = wdens0/(10.*sigmaw)
2405	!c sigmaw =max(sigmaw,sigd_con)
2406	!c sigmaw =max(sigmaw,sigmad)
2407
2408	!c calcul de la difference de vitesse verticale poche - zone non perturbee
2409
2410	z= 0.
2411	dp_deltomg(1:klev)=0.
2412	omg(1:klev+1)=0.
2413
2414	omg(1) = 0.
2415	dp_deltomg(1) = -(gflCstar)/(sigmaw (1-sigmaw))
2416
2417	DO k=2,ktop
2418	dz = -(Ph(k)-Ph(k-1))/(rho(k-1)*rg)
2419	z = z+dz
2420	dp_deltomg(k)= dp_deltomg(1)
2421	omg(k)= dp_deltomg(1)*z
2422	ENDDO
2423
2424	dztop=-(Ptop-Ph(ktop))/(rho(ktop)*rg)
2425	ztop = z+dztop
2426	omgtop=dp_deltomg(1)*ztop
2427
2428
2429	!c Conversion de la vitesse verticale de m/s a Pa/s
2430
2431	omgtop = -rho(ktop)rgomgtop
2432	dp_deltomg(1) = omgtop/(ptop-ph(1))
2433
2434	DO k = 1,ktop
2435	omg(k) = - rho(k)rgomg(k)
2436	dp_deltomg(k) = dp_deltomg(1)
2437	ENDDO
2438
2439	!c raccordement lineaire de omg de ptop a pupper
2440
2441	IF (kupper .GT. ktop) THEN
2442	omg(kupper+1) = - Rg*amdwn(kupper+1)/sigmaw &
2443	& + Rg*amup(kupper+1)/(1.-sigmaw)
2444	dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(Ptop-Pupper)
2445	DO k=ktop+1,kupper
2446	dp_deltomg(k) = dp_deltomg(kupper)
2447	omg(k) = omgtop+(ph(k)-Ptop)*dp_deltomg(kupper)
2448	ENDDO
2449	ENDIF
2450
2451	!c Compute wake average vertical velocity omgbw
2452
2453	DO k = 1,klev+1
2454	omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k)
2455	ENDDO
2456
2457	!c and its vertical gradient dp_omgbw
2458
2459	DO k = 1,klev
2460	dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
2461	ENDDO
2462
2463
2464	!c Upstream coefficients for omgb velocity
2465	!c-- (alpha_up(k) is the coefficient of the value at level k)
2466	!c-- (1-alpha_up(k) is the coefficient of the value at level k-1)
2467
2468	DO k = 1,klev
2469	alpha_up(k) = 0.
2470	IF (omgb(k) .GT. 0.) alpha_up(k) = 1.
2471	ENDDO
2472
2473	!c Matrix expressing [The,deltatw] from [Th1,Th2]
2474
2475	RRe1 = 1.-sigmaw
2476	RRe2 = sigmaw
2477	RRd1 = -1.
2478	RRd2 = 1.
2479
2480	!c Get [Th1,Th2], dth and [q1,q2]
2481
2482	DO k = 1,kupper+1
2483	dth(k) = deltatw(k)/ppi(k)
2484	Th1(k) = the(k) - sigmaw *dth(k) ! undisturbed area
2485	Th2(k) = the(k) + (1.-sigmaw)*dth(k) ! wake
2486	q1(k) = qe(k) - sigmaw *deltaqw(k) ! undisturbed area
2487	q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
2488	ENDDO
2489
2490	D_Th1(1) = 0.
2491	D_Th2(1) = 0.
2492	D_dth(1) = 0.
2493	D_q1(1) = 0.
2494	D_q2(1) = 0.
2495	D_dq(1) = 0.
2496
2497	DO k = 2,kupper+1 ! loop on interfaces
2498	D_Th1(k) = Th1(k-1)-Th1(k)
2499	D_Th2(k) = Th2(k-1)-Th2(k)
2500	D_dth(k) = dth(k-1)-dth(k)
2501	D_q1(k) = q1(k-1)-q1(k)
2502	D_q2(k) = q2(k-1)-q2(k)
2503	D_dq(k) = deltaqw(k-1)-deltaqw(k)
2504	ENDDO
2505
2506	omgbdth(1) = 0.
2507	omgbdq(1) = 0.
2508
2509	DO k = 2,kupper+1 ! loop on interfaces
2510	omgbdth(k) = omgb(k)*( dth(k-1) - dth(k))
2511	omgbdq(k) = omgb(k)*(deltaqw(k-1) - deltaqw(k))
2512	ENDDO
2513
2514
2515	!c-----------------------------------------------------------------
2516	DO k=1,kupper-1
2517	!c-----------------------------------------------------------------
2518	!c
2519	!c Compute redistribution (advective) term
2520	!c
2521	d_deltatw(k) = &
2522	& dtimesub/(Ph(k)-Ph(k+1))*( &
2523	& RRd1omg(k )sigmaw *D_Th1(k) &
2524	& -RRd2omg(k+1)(1.-sigmaw)*D_Th2(k+1) &
2525	& -(1.-alpha_up(k))omgbdth(k) - alpha_up(k+1)omgbdth(k+1) &
2526	& )*ppi(k)
2527	!c print*,'d_deltatw=',d_deltatw(k)
2528	!c
2529	d_deltaqw(k) = &
2530	& dtimesub/(Ph(k)-Ph(k+1))*( &
2531	& RRd1omg(k )sigmaw *D_q1(k) &
2532	& -RRd2omg(k+1)(1.-sigmaw)*D_q2(k+1) &
2533	& -(1.-alpha_up(k))omgbdq(k) - alpha_up(k+1)omgbdq(k+1) &
2534	& )
2535	!c print*,'d_deltaqw=',d_deltaqw(k)
2536	!c
2537	!c and increment large scale tendencies
2538	!c
2539	dtls(k) = dtls(k) + &
2540	& dtimesub*( &
2541	& ( RRe1omg(k )sigmaw *D_Th1(k) &
2542	& -RRe2omg(k+1)(1.-sigmaw)*D_Th2(k+1) ) &
2543	& /(Ph(k)-Ph(k+1)) &
2544	& -sigmaw(1.-sigmaw)dth(k)*dp_deltomg(k) &
2545	& )*ppi(k)
2546	!c print*,'dtls=',dtls(k)
2547	!c
2548	dqls(k) = dqls(k) + &
2549	& dtimesub*( &
2550	& ( RRe1omg(k )sigmaw *D_q1(k) &
2551	& -RRe2omg(k+1)(1.-sigmaw)*D_q2(k+1) ) &
2552	& /(Ph(k)-Ph(k+1)) &
2553	& -sigmaw(1.-sigmaw)deltaqw(k)*dp_deltomg(k) &
2554	& )
2555	!c print*,'dqls=',dqls(k)
2556
2557	!c-------------------------------------------------------------------
2558	ENDDO
2559	!c------------------------------------------------------------------
2560
2561	!C Increment state variables
2562
2563	DO k = 1,kupper-1
2564
2565	!c Coefficient de répartition
2566
2567	Crep(k)=Crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
2568	Crep(k)=Crep(k)+Crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
2569
2570
2571	!c Reintroduce compensating subsidence term.
2572
2573	!c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
2574	!c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
2575	!c . /(1-sigmaw)
2576	!c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
2577	!c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
2578	!c . /(1-sigmaw)
2579	!c
2580	!c dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
2581	!c dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
2582	!c . /(1-sigmaw)
2583	!c dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
2584	!c dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
2585	!c . /(1-sigmaw)
2586
2587	dtKE(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw))
2588	dqKE(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw))
2589	!c print*,'dtKE=',dtKE(k)
2590	!c print*,'dqKE=',dqKE(k)
2591	!c
2592	dtPBL(k)=(wdtPBL(k)/sigmaw - udtPBL(k)/(1.-sigmaw))
2593	dqPBL(k)=(wdqPBL(k)/sigmaw - udqPBL(k)/(1.-sigmaw))
2594	!c
2595	spread(k) = (1.-sigmaw)dp_deltomg(k)+gflCstar/sigmaw
2596	!c print*,'spread=',spread(k)
2597
2598
2599	!c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
2600
2601	d_deltat_gw(k)=d_deltat_gw(k)-Tgw(k)deltatw(k) dtimesub
2602	!c print*,'d_delta_gw=',d_deltat_gw(k)
2603	ff=d_deltatw(k)/dtimesub
2604
2605	!c Sans GW
2606	!c
2607	!c deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
2608	!c
2609	!c GW formule 1
2610	!c
2611	!c deltatw(k) = deltatw(k)+dtimesub*
2612	!c $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
2613	!c
2614	!c GW formule 2
2615
2616	IF (dtimesub*Tgw(k).lt.1.e-10) THEN
2617	deltatw(k) = deltatw(k)+dtimesub* &
2618	& (ff+dtKE(k)+dtPBL(k) &
2619	& - spread(k)deltatw(k)-Tgw(k)deltatw(k))
2620	ELSE
2621	deltatw(k) = deltatw(k)+1/Tgw(k)(1-exp(-dtimesubTgw(k)))* &
2622	& (ff+dtKE(k)+dtPBL(k) &
2623	& - spread(k)deltatw(k)-Tgw(k)deltatw(k))
2624	ENDIF
2625
2626	dth(k) = deltatw(k)/ppi(k)
2627
2628	gg=d_deltaqw(k)/dtimesub
2629
2630	deltaqw(k) = deltaqw(k) + &
2631	& dtimesub(gg+ dqKE(k)+dqPBL(k) - spread(k)deltaqw(k))
2632
2633	d_deltatw2(k)=d_deltatw2(k)+d_deltatw(k)
2634	d_deltaqw2(k)=d_deltaqw2(k)+d_deltaqw(k)
2635	ENDDO
2636
2637	!C And update large scale variables
2638
2639	DO k = 1,kupper
2640	te(k) = te0(k) + dtls(k)
2641	qe(k) = qe0(k) + dqls(k)
2642	the(k) = te(k)/ppi(k)
2643	ENDDO
2644
2645	!c Determine Ptop from buoyancy integral
2646	!c----------------------------------------------------------------------
2647
2648	!c-1/ Pressure of the level where dth changes sign.
2649
2650	Ptop_provis=ph(1)
2651
2652	DO k= 2,klev
2653	IF (dth(k) .GT. -delta_t_min .and. &
2654	& dth(k-1).LT. -delta_t_min) THEN
2655	Ptop_provis = ((dth(k)+delta_t_min)*p(k-1) &
2656	& - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2657	GO TO 65
2658	ENDIF
2659	ENDDO
2660	65 CONTINUE
2661
2662	!c-2/ dth integral
2663
2664	sum_dth = 0.
2665	dthmin = -delta_t_min
2666	z = 0.
2667
2668	DO k = 1,klev
2669	dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
2670	IF (dz .le. 0) GO TO 70
2671	z = z+dz
2672	sum_dth = sum_dth + dth(k)*dz
2673	dthmin = min(dthmin,dth(k))
2674	ENDDO
2675	70 CONTINUE
2676
2677	!c-3/ height of triangle with area= sum_dth and base = dthmin
2678
2679	hw = 2.*sum_dth/min(dthmin,-0.5)
2680	hw = max(hwmin,hw)
2681
2682	!c-4/ now, get Ptop
2683
2684	ktop = 0
2685	z=0.
2686
2687	DO k = 1,klev
2688	dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw-z)
2689	IF (dz .le. 0) GO TO 75
2690	z = z+dz
2691	Ptop = Ph(k)-rho(k)rgdz
2692	ktop = k
2693	ENDDO
2694	75 CONTINUE
2695
2696	!c-5/Correct ktop and ptop
2697
2698	Ptop_new=ptop
2699
2700	DO k= ktop,2,-1
2701	IF (dth(k) .GT. -delta_t_min .and. &
2702	& dth(k-1).LT. -delta_t_min) THEN
2703	Ptop_new = ((dth(k)+delta_t_min)*p(k-1) &
2704	& - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2705	GO TO 275
2706	ENDIF
2707	ENDDO
2708	275 CONTINUE
2709
2710	ptop = ptop_new
2711
2712	DO k=klev,1,-1
2713	IF (ph(k+1) .LT. ptop) ktop=k
2714	ENDDO
2715
2716	!c-6/ Set deltatw & deltaqw to 0 above kupper
2717
2718	DO k = kupper,klev
2719	deltatw(k) = 0.
2720	deltaqw(k) = 0.
2721	ENDDO
2722
2723	!c------------------------------------------------------------------
2724	ENDDO ! end sub-timestep loop
2725	!C -----------------------------------------------------------------
2726
2727	!c Get back to tendencies per second
2728
2729	DO k = 1,kupper-1
2730	dtls(k) = dtls(k)/dtime
2731	dqls(k) = dqls(k)/dtime
2732	d_deltatw2(k)=d_deltatw2(k)/dtime
2733	d_deltaqw2(k)=d_deltaqw2(k)/dtime
2734	d_deltat_gw(k) = d_deltat_gw(k)/dtime
2735	ENDDO
2736
2737	!C 2.1 - Undisturbed area and Wake integrals
2738	!C ---------------------------------------------------------
2739
2740	z = 0.
2741	sum_thu = 0.
2742	sum_tu = 0.
2743	sum_qu = 0.
2744	sum_thvu = 0.
2745	sum_dth = 0.
2746	sum_dq = 0.
2747	sum_rho = 0.
2748	sum_dtdwn = 0.
2749	sum_dqdwn = 0.
2750
2751	av_thu = 0.
2752	av_tu =0.
2753	av_qu =0.
2754	av_thvu = 0.
2755	av_dth = 0.
2756	av_dq = 0.
2757	av_rho =0.
2758	av_dtdwn =0.
2759	av_dqdwn = 0.
2760
2761	!C Potential temperatures and humidity
2762	!c----------------------------------------------------------
2763
2764	DO k =1,klev
2765	rho(k) = p(k)/(rd*te(k))
2766	IF(k .eq. 1) THEN
2767	rhoh(k) = ph(k)/(rd*te(k))
2768	zhh(k)=0
2769	ELSE
2770	rhoh(k) = ph(k)2./(rd(te(k)+te(k-1)))
2771	zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
2772	ENDIF
2773	the(k) = te(k)/ppi(k)
2774	thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
2775	tu(k) = te(k) - deltatw(k)*sigmaw
2776	qu(k) = qe(k) - deltaqw(k)*sigmaw
2777	rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
2778	dth(k) = deltatw(k)/ppi(k)
2779
2780	ENDDO
2781
2782	!C Integrals (and wake top level number)
2783	!C -----------------------------------------------------------
2784
2785	!C Initialize sum_thvu to 1st level virt. pot. temp.
2786
2787	z = 1.
2788	dz = 1.
2789	sum_thvu = thu(1)(1.+epsqu(1))*dz
2790	sum_dth = 0.
2791
2792	DO k = 1,klev
2793	dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
2794
2795	IF (dz .LE. 0) GO TO 51
2796	z = z+dz
2797	sum_thu = sum_thu + thu(k)*dz
2798	sum_tu = sum_tu + tu(k)*dz
2799	sum_qu = sum_qu + qu(k)*dz
2800	sum_thvu = sum_thvu + thu(k)(1.+epsqu(k))*dz
2801	sum_dth = sum_dth + dth(k)*dz
2802	sum_dq = sum_dq + deltaqw(k)*dz
2803	sum_rho = sum_rho + rhow(k)*dz
2804	sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
2805	sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
2806	ENDDO
2807	51 CONTINUE
2808
2809	hw0 = z
2810
2811	!C 2.1 - WAPE and mean forcing computation
2812	!C-------------------------------------------------------------
2813
2814	!C Means
2815
2816	av_thu = sum_thu/hw0
2817	av_tu = sum_tu/hw0
2818	av_qu = sum_qu/hw0
2819	av_thvu = sum_thvu/hw0
2820	av_dth = sum_dth/hw0
2821	av_dq = sum_dq/hw0
2822	av_rho = sum_rho/hw0
2823	av_dtdwn = sum_dtdwn/hw0
2824	av_dqdwn = sum_dqdwn/hw0
2825
2826	wape2 = - rghw0(av_dth &
2827	& + eps(av_thuav_dq+av_dthav_qu+av_dthav_dq ))/av_thvu
2828
2829
2830	!C 2.2 Prognostic variable update
2831	!C ------------------------------------------------------------
2832
2833	!C Filter out bad wakes
2834
2835	IF ( wape2 .LT. 0.) THEN
2836	if(prt_level.ge.10) print*,'wape2<0'
2837	wape2 = 0.
2838	hw = hwmin
2839	sigmaw = max(sigmad,sigd_con)
2840	fip = 0.
2841	DO k = 1,klev
2842	deltatw(k) = 0.
2843	deltaqw(k) = 0.
2844	dth(k) = 0.
2845	ENDDO
2846	ELSE
2847	if(prt_level.ge.10) print*,'wape2>0'
2848	Cstar2 = starksqrt(2.wape2)
2849
2850	ENDIF
2851
2852	ktopw = ktop
2853
2854	IF (ktopw .gt. 0) then
2855
2856	!Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
2857	!ccc heff = 600.
2858	!C Utilisation de la hauteur hw
2859	!cc heff = 0.7*hw
2860	heff = hw
2861
2862	FIP = 0.5rho(ktopw)Cstar2*3heff2sqrt(sigmawwdens3.14)
2863	FIP = alpk * FIP
2864	!Cjyg2
2865	ELSE
2866	FIP = 0.
2867	ENDIF
2868
2869
2870	!C Limitation de sigmaw
2871	!c
2872	!C sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
2873	!C alors il disparait en se mélangeant à la partie undisturbed
2874
2875	! correction NICOLAS . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
2876	IF ((sigmaw.GT.0.9).or. &
2877	& ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN
2878	!IM cf NR/JYG 251108 . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
2879	!c IF (sigmaw.GT.0.9) THEN
2880	DO k = 1,klev
2881	dtls(k) = 0.
2882	dqls(k) = 0.
2883	deltatw(k) = 0.
2884	deltaqw(k) = 0.
2885	ENDDO
2886	wape = 0.
2887	hw = hwmin
2888	sigmaw = sigmad
2889	fip = 0.
2890	ENDIF
2891
2892	RETURN
2893	END SUBROUTINE WAKE_scal
2894
2895
2896

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

Download in other formats:

Original Format