Context Navigation

wake.F @ 1233

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

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

Download in other formats:

Original Format