Context Navigation

wake.F @ 3400

Last change on this file since 3400 was 1277, checked in by musat, 15 years ago
Corrections conservation de l'eau dans les descentes insaturees pour la nouvelle physique JYG/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 77.3 KB

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

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

Download in other formats:

Original Format