Context Navigation

wake.F @ 1127

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

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

Download in other formats:

Original Format