Context Navigation

wake.F @ 357

Last change on this file since 357 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 82.8 KB

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

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

Download in other formats:

Original Format