Context Navigation

wake.F @ 5456

Last change on this file since 5456 was 1861, checked in by lguez, 11 years ago
Replaced #include by include. (Use as little C preprocessing as possible.)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 82.8 KB

Line
1	!
2	! $Id: wake.F 1861 2013-09-10 09:55:26Z aborella $
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	o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
9	o ,dtls,dqls
10	o ,ktopw,omgbdth,dp_omgb,wdens
11	o ,tu,qu
12	o ,dtKE,dqKE
13	o ,dtPBL,dqPBL
14	o ,omg,dp_deltomg,spread
15	o ,Cstar,d_deltat_gw
16	o ,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	cIM 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	cIM 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	cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
803	cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
804	cIM 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	cIM 060208 manque DO i + remplace DO k=1,kupper(i)
1209	cIM 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	cIM 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
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
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	o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
1837	o ,dtls,dqls
1838	o ,ktopw,omgbdth,dp_omgb,wdens
1839	o ,tu,qu
1840	o ,dtKE,dqKE
1841	o ,dtPBL,dqPBL
1842	o ,omg,dp_deltomg,spread
1843	o ,Cstar,d_deltat_gw
1844	o ,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	cIM 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
2894
2895
2896

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

Download in other formats:

Original Format