Context Navigation

wake.F90 @ 2641

Last change on this file since 2641 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.9 KB

Line
1
2	! $Id: wake.F90 2641 2016-09-29 21:26:46Z fairhead $
3
4	SUBROUTINE wake(p, ph, pi, dtime, &
5	te0, qe0, omgb, &
6	dtdwn, dqdwn, amdwn, amup, dta, dqa, &
7	sigd_con, &
8	deltatw, deltaqw, sigmaw, wdens, & ! state variables
9	dth, hw, wape, fip, gfl, &
10	dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
11	dtke, dqke, omg, dp_deltomg, spread, cstar, &
12	d_deltat_gw, &
13	d_deltatw2, d_deltaqw2, d_sigmaw2, d_wdens2) ! tendencies
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
25	USE ioipsl_getin_p_mod, ONLY : getin_p
26	USE dimphy
27	use mod_phys_lmdz_para
28	USE print_control_mod, ONLY: prt_level
29	IMPLICIT NONE
30	! ============================================================================
31
32
33	! But : Decrire le comportement des poches froides apparaissant dans les
34	! grands systemes convectifs, et fournir l'energie disponible pour
35	! le declenchement de nouvelles colonnes convectives.
36
37	! State variables :
38	! deltatw : temperature difference between wake and off-wake regions
39	! deltaqw : specific humidity difference between wake and off-wake regions
40	! sigmaw : fractional area covered by wakes.
41	! wdens : number of wakes per unit area
42
43	! Variable de sortie :
44
45	! wape : WAke Potential Energy
46	! fip : Front Incident Power (W/m2) - ALP
47	! gfl : Gust Front Length per unit area (m-1)
48	! dtls : large scale temperature tendency due to wake
49	! dqls : large scale humidity tendency due to wake
50	! hw : hauteur de la poche
51	! dp_omgb : vertical gradient of large scale omega
52	! wdens : densite de poches
53	! omgbdth: flux of Delta_Theta transported by LS omega
54	! dtKE : differential heating (wake - unpertubed)
55	! dqKE : differential moistening (wake - unpertubed)
56	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
57	! dp_deltomg : vertical gradient of omg (s-1)
58	! spread : spreading term in d_t_wake and d_q_wake
59	! deltatw : updated temperature difference (T_w-T_u).
60	! deltaqw : updated humidity difference (q_w-q_u).
61	! sigmaw : updated wake fractional area.
62	! d_deltat_gw : delta T tendency due to GW
63
64	! Variables d'entree :
65
66	! aire : aire de la maille
67	! te0 : temperature dans l'environnement (K)
68	! qe0 : humidite dans l'environnement (kg/kg)
69	! omgb : vitesse verticale moyenne sur la maille (Pa/s)
70	! dtdwn: source de chaleur due aux descentes (K/s)
71	! dqdwn: source d'humidite due aux descentes (kg/kg/s)
72	! dta : source de chaleur due courants satures et detrain (K/s)
73	! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
74	! amdwn: flux de masse total des descentes, par unite de
75	! surface de la maille (kg/m2/s)
76	! amup : flux de masse total des ascendances, par unite de
77	! surface de la maille (kg/m2/s)
78	! p : pressions aux milieux des couches (Pa)
79	! ph : pressions aux interfaces (Pa)
80	! pi : (p/p_0)**kapa (adim)
81	! dtime: increment temporel (s)
82
83	! Variables internes :
84
85	! rhow : masse volumique de la poche froide
86	! rho : environment density at P levels
87	! rhoh : environment density at Ph levels
88	! te : environment temperature \| may change within
89	! qe : environment humidity \| sub-time-stepping
90	! the : environment potential temperature
91	! thu : potential temperature in undisturbed area
92	! tu : temperature in undisturbed area
93	! qu : humidity in undisturbed area
94	! dp_omgb: vertical gradient og LS omega
95	! omgbw : wake average vertical omega
96	! dp_omgbw: vertical gradient of omgbw
97	! omgbdq : flux of Delta_q transported by LS omega
98	! dth : potential temperature diff. wake-undist.
99	! th1 : first pot. temp. for vertical advection (=thu)
100	! th2 : second pot. temp. for vertical advection (=thw)
101	! q1 : first humidity for vertical advection
102	! q2 : second humidity for vertical advection
103	! d_deltatw : terme de redistribution pour deltatw
104	! d_deltaqw : terme de redistribution pour deltaqw
105	! deltatw0 : deltatw initial
106	! deltaqw0 : deltaqw initial
107	! hw0 : hw initial
108	! sigmaw0: sigmaw initial
109	! amflux : horizontal mass flux through wake boundary
110	! wdens_ref: initial number of wakes per unit area (3D) or per
111	! unit length (2D), at the beginning of each time step
112	! Tgw : 1 sur la période de onde de gravité
113	! Cgw : vitesse de propagation de onde de gravité
114	! LL : distance entre 2 poches
115
116	! -------------------------------------------------------------------------
117	! Déclaration de variables
118	! -------------------------------------------------------------------------
119
120	include "YOMCST.h"
121	include "cvthermo.h"
122
123	! Arguments en entree
124	! --------------------
125
126	REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi
127	REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph, omgb
128	REAL, INTENT(IN) :: dtime
129	REAL, DIMENSION (klon, klev), INTENT(IN) :: te0, qe0
130	REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn
131	REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup
132	REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa
133	REAL, DIMENSION (klon), INTENT(IN) :: sigd_con
134
135	!
136	! Input/Output
137	! State variables
138	REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw
139	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw
140	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens
141
142	! Sorties
143	! --------
144
145	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth
146	REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu
147	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls
148	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke
149	REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread
150	REAL, DIMENSION (klon, klev+1), INTENT(OUT) :: omgbdth, omg
151	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg
152	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw
153	REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar
154	INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw
155	! Tendencies of state variables
156	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2
157	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_wdens2
158
159	! Variables internes
160	! -------------------
161
162	! Variables à fixer
163	REAL :: alon
164	LOGICAL, SAVE :: first = .TRUE.
165	!$OMP THREADPRIVATE(first)
166	REAL, SAVE :: stark, wdens_ref, coefgw, alpk
167	REAL, SAVE :: crep_upper, crep_sol
168	!$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
169
170	REAL :: delta_t_min
171	INTEGER :: nsub
172	REAL :: dtimesub
173	REAL :: sigmad, hwmin, wapecut
174	REAL :: sigmaw_max
175	REAL :: dens_rate
176	REAL :: wdens0
177	! IM 080208
178	LOGICAL, DIMENSION (klon) :: gwake
179
180	! Variables de sauvegarde
181	REAL, DIMENSION (klon, klev) :: deltatw0
182	REAL, DIMENSION (klon, klev) :: deltaqw0
183	REAL, DIMENSION (klon, klev) :: te, qe
184	REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
185
186	! Variables pour les GW
187	REAL, DIMENSION (klon) :: ll
188	REAL, DIMENSION (klon, klev) :: n2
189	REAL, DIMENSION (klon, klev) :: cgw
190	REAL, DIMENSION (klon, klev) :: tgw
191
192	! Variables liées au calcul de hw
193	REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
194	REAL, DIMENSION (klon) :: sum_dth
195	REAL, DIMENSION (klon) :: dthmin
196	REAL, DIMENSION (klon) :: z, dz, hw0
197	INTEGER, DIMENSION (klon) :: ktop, kupper
198
199	! Sub-timestep tendencies and related variables
200	REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw
201	REAL, DIMENSION (klon, klev) :: d_te, d_qe
202	REAL, DIMENSION (klon) :: d_sigmaw, alpha
203	REAL, DIMENSION (klon) :: q0_min, q1_min
204	LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw
205	REAL, SAVE :: epsilon
206	!$OMP THREADPRIVATE(epsilon)
207	DATA epsilon/1.E-15/
208
209	! Autres variables internes
210	INTEGER ::isubstep, k, i
211
212	REAL :: sigmaw_targ
213
214	REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
215	REAL, DIMENSION (klon) :: sum_dq, sum_rho
216	REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
217	REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
218	REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
219	REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
220
221	REAL, DIMENSION (klon, klev) :: rho, rhow
222	REAL, DIMENSION (klon, klev+1) :: rhoh
223	REAL, DIMENSION (klon, klev) :: rhow_moyen
224	REAL, DIMENSION (klon, klev) :: zh
225	REAL, DIMENSION (klon, klev+1) :: zhh
226	REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
227
228	REAL, DIMENSION (klon, klev) :: the, thu
229
230	REAL, DIMENSION (klon, klev+1) :: omgbw
231	REAL, DIMENSION (klon) :: pupper
232	REAL, DIMENSION (klon) :: omgtop
233	REAL, DIMENSION (klon, klev) :: dp_omgbw
234	REAL, DIMENSION (klon) :: ztop, dztop
235	REAL, DIMENSION (klon, klev) :: alpha_up
236
237	REAL, DIMENSION (klon) :: rre1, rre2
238	REAL :: rrd1, rrd2
239	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
240	REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
241	REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
242	REAL, DIMENSION (klon, klev) :: omgbdq
243
244	REAL, DIMENSION (klon) :: ff, gg
245	REAL, DIMENSION (klon) :: wape2, cstar2, heff
246
247	REAL, DIMENSION (klon, klev) :: crep
248
249	REAL, DIMENSION (klon, klev) :: ppi
250
251	! cc nrlmd
252	REAL, DIMENSION (klon) :: death_rate, nat_rate
253	REAL, DIMENSION (klon, klev) :: entr
254	REAL, DIMENSION (klon, klev) :: detr
255
256	REAL, DIMENSION(klon) :: sigmaw_in ! pour les prints
257
258	! -------------------------------------------------------------------------
259	! Initialisations
260	! -------------------------------------------------------------------------
261
262	! print*, 'wake initialisations'
263
264	! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
265	! -------------------------------------------------------------------------
266
267	DATA wapecut, sigmad, hwmin/5., .02, 10./
268	! cc nrlmd
269	DATA sigmaw_max/0.4/
270	DATA dens_rate/0.1/
271	! cc
272	! Longueur de maille (en m)
273	! -------------------------------------------------------------------------
274
275	! ALON = 3.e5
276	alon = 1.E6
277
278
279	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
280
281	! coefgw : Coefficient pour les ondes de gravité
282	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
283	! wdens : Densité de poche froide par maille
284	! -------------------------------------------------------------------------
285
286	! cc nrlmd coefgw=10
287	! coefgw=1
288	! wdens0 = 1.0/(alon**2)
289	! cc nrlmd wdens = 1.0/(alon**2)
290	! cc nrlmd stark = 0.50
291	! CRtest
292	! cc nrlmd alpk=0.1
293	! alpk = 1.0
294	! alpk = 0.5
295	! alpk = 0.05
296
297	if (first) then
298	crep_upper = 0.9
299	crep_sol = 1.0
300
301	! cc nrlmd Lecture du fichier wake_param.data
302	stark=0.33
303	CALL getin_p('stark',stark)
304	alpk=0.25
305	CALL getin_p('alpk',alpk)
306	wdens_ref=8.E-12
307	CALL getin_p('wdens_ref',wdens_ref)
308	coefgw=4.
309	CALL getin_p('coefgw',coefgw)
310
311	WRITE(,) 'stark=', stark
312	WRITE(,) 'alpk=', alpk
313	WRITE(,) 'wdens_ref=', wdens_ref
314	WRITE(,) 'coefgw=', coefgw
315
316	first=.false.
317	endif
318
319	! Initialisation de toutes des densites a wdens_ref.
320	! Les densites peuvent evoluer si les poches debordent
321	! (voir au tout debut de la boucle sur les substeps)
322	wdens(:) = wdens_ref
323
324	! print*,'stark',stark
325	! print*,'alpk',alpk
326	! print*,'wdens',wdens
327	! print*,'coefgw',coefgw
328	! cc
329	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
330	! -------------------------------------------------------------------------
331
332	delta_t_min = 0.2
333
334	! 1. - Save initial values and initialize tendencies
335	! --------------------------------------------------
336
337	DO k = 1, klev
338	DO i = 1, klon
339	ppi(i, k) = pi(i, k)
340	deltatw0(i, k) = deltatw(i, k)
341	deltaqw0(i, k) = deltaqw(i, k)
342	te(i, k) = te0(i, k)
343	qe(i, k) = qe0(i, k)
344	dtls(i, k) = 0.
345	dqls(i, k) = 0.
346	d_deltat_gw(i, k) = 0.
347	d_te(i, k) = 0.
348	d_qe(i, k) = 0.
349	d_deltatw(i, k) = 0.
350	d_deltaqw(i, k) = 0.
351	! IM 060508 beg
352	d_deltatw2(i, k) = 0.
353	d_deltaqw2(i, k) = 0.
354	! IM 060508 end
355	END DO
356	END DO
357	DO i = 1, klon
358	sigmaw_in(i) = sigmaw(i)
359	END DO
360	! sigmaw1=sigmaw
361	! IF (sigd_con.GT.sigmaw1) THEN
362	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
363	! ENDIF
364	DO i = 1, klon
365	! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
366	!jyg<
367	!! sigmaw(i) = amax1(sigmaw(i), sigmad)
368	!! sigmaw(i) = amin1(sigmaw(i), 0.99)
369	sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
370	d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
371	sigmaw(i) = sigmaw_targ
372	!>jyg
373	sigmaw0(i) = sigmaw(i)
374	wape(i) = 0.
375	wape2(i) = 0.
376	d_sigmaw(i) = 0.
377	d_wdens2(i) = 0.
378	ktopw(i) = 0
379	END DO
380
381
382	! 2. - Prognostic part
383	! --------------------
384
385
386	! 2.1 - Undisturbed area and Wake integrals
387	! ---------------------------------------------------------
388
389	DO i = 1, klon
390	z(i) = 0.
391	ktop(i) = 0
392	kupper(i) = 0
393	sum_thu(i) = 0.
394	sum_tu(i) = 0.
395	sum_qu(i) = 0.
396	sum_thvu(i) = 0.
397	sum_dth(i) = 0.
398	sum_dq(i) = 0.
399	sum_rho(i) = 0.
400	sum_dtdwn(i) = 0.
401	sum_dqdwn(i) = 0.
402
403	av_thu(i) = 0.
404	av_tu(i) = 0.
405	av_qu(i) = 0.
406	av_thvu(i) = 0.
407	av_dth(i) = 0.
408	av_dq(i) = 0.
409	av_rho(i) = 0.
410	av_dtdwn(i) = 0.
411	av_dqdwn(i) = 0.
412	END DO
413
414	! Distance between wakes
415	DO i = 1, klon
416	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
417	END DO
418	! Potential temperatures and humidity
419	! ----------------------------------------------------------
420	DO k = 1, klev
421	DO i = 1, klon
422	! write(,)'wake 1',i,k,rd,te(i,k)
423	rho(i, k) = p(i, k)/(rd*te(i,k))
424	! write(,)'wake 2',rho(i,k)
425	IF (k==1) THEN
426	! write(,)'wake 3',i,k,rd,te(i,k)
427	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
428	! write(,)'wake 4',i,k,rd,te(i,k)
429	zhh(i, k) = 0
430	ELSE
431	! write(,)'wake 5',rd,(te(i,k)+te(i,k-1))
432	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
433	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
434	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
435	END IF
436	! write(,)'wake 7',ppi(i,k)
437	the(i, k) = te(i, k)/ppi(i, k)
438	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
439	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
440	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
441	! write(,)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
442	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
443	dth(i, k) = deltatw(i, k)/ppi(i, k)
444	END DO
445	END DO
446
447	DO k = 1, klev - 1
448	DO i = 1, klon
449	IF (k==1) THEN
450	n2(i, k) = 0
451	ELSE
452	n2(i, k) = amax1(0., -rg*2/the(i,k)rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
453	(p(i,k+1)-p(i,k-1)))
454	END IF
455	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
456
457	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
458	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
459	END DO
460	END DO
461
462	DO i = 1, klon
463	n2(i, klev) = 0
464	zh(i, klev) = 0
465	cgw(i, klev) = 0
466	tgw(i, klev) = 0
467	END DO
468
469	! Calcul de la masse volumique moyenne de la colonne (bdlmd)
470	! -----------------------------------------------------------------
471
472	DO k = 1, klev
473	DO i = 1, klon
474	epaisseur1(i, k) = 0.
475	epaisseur2(i, k) = 0.
476	END DO
477	END DO
478
479	DO i = 1, klon
480	epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
481	epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
482	rhow_moyen(i, 1) = rhow(i, 1)
483	END DO
484
485	DO k = 2, klev
486	DO i = 1, klon
487	epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
488	epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
489	rhow_moyen(i, k) = (rhow_moyen(i,k-1)epaisseur2(i,k-1)+rhow(i,k) &
490	epaisseur1(i,k))/epaisseur2(i, k)
491	END DO
492	END DO
493
494
495	! Choose an integration bound well above wake top
496	! -----------------------------------------------------------------
497
498	! Pupper = 50000. ! melting level
499	! Pupper = 60000.
500	! Pupper = 80000. ! essais pour case_e
501	DO i = 1, klon
502	pupper(i) = 0.6*ph(i, 1)
503	pupper(i) = max(pupper(i), 45000.)
504	! cc Pupper(i) = 60000.
505	END DO
506
507
508	! Determine Wake top pressure (Ptop) from buoyancy integral
509	! --------------------------------------------------------
510
511	! -1/ Pressure of the level where dth becomes less than delta_t_min.
512
513	DO i = 1, klon
514	ptop_provis(i) = ph(i, 1)
515	END DO
516	DO k = 2, klev
517	DO i = 1, klon
518
519	! IM v3JYG; ptop_provis(i).LT. ph(i,1)
520
521	IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
522	ptop_provis(i)==ph(i,1)) THEN
523	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &
524	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
525	END IF
526	END DO
527	END DO
528
529	! -2/ dth integral
530
531	DO i = 1, klon
532	sum_dth(i) = 0.
533	dthmin(i) = -delta_t_min
534	z(i) = 0.
535	END DO
536
537	DO k = 1, klev
538	DO i = 1, klon
539	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
540	IF (dz(i)>0) THEN
541	z(i) = z(i) + dz(i)
542	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
543	dthmin(i) = amin1(dthmin(i), dth(i,k))
544	END IF
545	END DO
546	END DO
547
548	! -3/ height of triangle with area= sum_dth and base = dthmin
549
550	DO i = 1, klon
551	hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
552	hw0(i) = amax1(hwmin, hw0(i))
553	END DO
554
555	! -4/ now, get Ptop
556
557	DO i = 1, klon
558	z(i) = 0.
559	ptop(i) = ph(i, 1)
560	END DO
561
562	DO k = 1, klev
563	DO i = 1, klon
564	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
565	IF (dz(i)>0) THEN
566	z(i) = z(i) + dz(i)
567	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
568	END IF
569	END DO
570	END DO
571
572
573	! -5/ Determination de ktop et kupper
574
575	DO k = klev, 1, -1
576	DO i = 1, klon
577	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
578	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
579	END DO
580	END DO
581
582	! On evite kupper = 1 et kupper = klev
583	DO i = 1, klon
584	kupper(i) = max(kupper(i), 2)
585	kupper(i) = min(kupper(i), klev-1)
586	END DO
587
588
589	! -6/ Correct ktop and ptop
590
591	DO i = 1, klon
592	ptop_new(i) = ptop(i)
593	END DO
594	DO k = klev, 2, -1
595	DO i = 1, klon
596	IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
597	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
598	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
599	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
600	END IF
601	END DO
602	END DO
603
604	DO i = 1, klon
605	ptop(i) = ptop_new(i)
606	END DO
607
608	DO k = klev, 1, -1
609	DO i = 1, klon
610	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
611	END DO
612	END DO
613
614	! -5/ Set deltatw & deltaqw to 0 above kupper
615
616	DO k = 1, klev
617	DO i = 1, klon
618	IF (k>=kupper(i)) THEN
619	deltatw(i, k) = 0.
620	deltaqw(i, k) = 0.
621	d_deltatw2(i,k) = -deltatw0(i,k)
622	d_deltaqw2(i,k) = -deltaqw0(i,k)
623	END IF
624	END DO
625	END DO
626
627
628	! Vertical gradient of LS omega
629
630	DO k = 1, klev
631	DO i = 1, klon
632	IF (k<=kupper(i)) THEN
633	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
634	END IF
635	END DO
636	END DO
637
638	! Integrals (and wake top level number)
639	! --------------------------------------
640
641	! Initialize sum_thvu to 1st level virt. pot. temp.
642
643	DO i = 1, klon
644	z(i) = 1.
645	dz(i) = 1.
646	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
647	sum_dth(i) = 0.
648	END DO
649
650	DO k = 1, klev
651	DO i = 1, klon
652	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
653	IF (dz(i)>0) THEN
654	z(i) = z(i) + dz(i)
655	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
656	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
657	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
658	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
659	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
660	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
661	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
662	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
663	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
664	END IF
665	END DO
666	END DO
667
668	DO i = 1, klon
669	hw0(i) = z(i)
670	END DO
671
672
673	! 2.1 - WAPE and mean forcing computation
674	! ---------------------------------------
675
676	! ---------------------------------------
677
678	! Means
679
680	DO i = 1, klon
681	av_thu(i) = sum_thu(i)/hw0(i)
682	av_tu(i) = sum_tu(i)/hw0(i)
683	av_qu(i) = sum_qu(i)/hw0(i)
684	av_thvu(i) = sum_thvu(i)/hw0(i)
685	! av_thve = sum_thve/hw0
686	av_dth(i) = sum_dth(i)/hw0(i)
687	av_dq(i) = sum_dq(i)/hw0(i)
688	av_rho(i) = sum_rho(i)/hw0(i)
689	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
690	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
691
692	wape(i) = -rghw0(i)(av_dth(i)+ &
693	epsim1(av_thu(i)av_dq(i)+av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
694
695	END DO
696
697	! 2.2 Prognostic variable update
698	! ------------------------------
699
700	! Filter out bad wakes
701
702	DO k = 1, klev
703	DO i = 1, klon
704	IF (wape(i)<0.) THEN
705	deltatw(i, k) = 0.
706	deltaqw(i, k) = 0.
707	dth(i, k) = 0.
708	d_deltatw2(i,k) = -deltatw0(i,k)
709	d_deltaqw2(i,k) = -deltaqw0(i,k)
710	END IF
711	END DO
712	END DO
713
714	DO i = 1, klon
715	IF (wape(i)<0.) THEN
716	wape(i) = 0.
717	cstar(i) = 0.
718	hw(i) = hwmin
719	!jyg<
720	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
721	sigmaw_targ = max(sigmad, sigd_con(i))
722	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
723	sigmaw(i) = sigmaw_targ
724	!>jyg
725	fip(i) = 0.
726	gwake(i) = .FALSE.
727	ELSE
728	cstar(i) = starksqrt(2.wape(i))
729	gwake(i) = .TRUE.
730	END IF
731	END DO
732
733
734	! Check qx and qw positivity
735	! --------------------------
736	DO i = 1, klon
737	q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), &
738	(qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
739	END DO
740	DO k = 2, klev
741	DO i = 1, klon
742	q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
743	(qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
744	IF (q1_min(i)<=q0_min(i)) THEN
745	q0_min(i) = q1_min(i)
746	END IF
747	END DO
748	END DO
749
750	DO i = 1, klon
751	ok_qx_qw(i) = q0_min(i) >= 0.
752	alpha(i) = 1.
753	END DO
754
755	! C -----------------------------------------------------------------
756	! Sub-time-stepping
757	! -----------------
758
759	nsub = 10
760	dtimesub = dtime/nsub
761
762	! ------------------------------------------------------------
763	DO isubstep = 1, nsub
764	! ------------------------------------------------------------
765
766	! wk_adv is the logical flag enabling wake evolution in the time advance
767	! loop
768	DO i = 1, klon
769	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
770	END DO
771
772	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
773	! négatif de ktop à kupper --------
774	! cc On calcule pour cela une densité wdens0 pour laquelle on
775	! aurait un entrainement nul ---
776	DO i = 1, klon
777	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
778	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
779	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
780	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
781	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
782	wdens0 = (sigmaw(i)/(4.3.14)) &
783	((1.-sigmaw(i))omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))cstar(i)))**(2)
784	IF (wdens(i)<=wdens0*1.1) THEN
785	wdens(i) = wdens0
786	END IF
787	! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
788	! c $ ,ph(i,1)-pupper(i)',
789	! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
790	! c $ ,ph(i,1)-pupper(i)
791	END IF
792	END DO
793
794	! cc nrlmd
795
796	DO i = 1, klon
797	IF (wk_adv(i)) THEN
798	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
799	!jyg<
800	!! sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
801	sigmaw_targ = min(sigmaw(i), sigmaw_max)
802	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
803	sigmaw(i) = sigmaw_targ
804	!>jyg
805	END IF
806	END DO
807
808	DO i = 1, klon
809	IF (wk_adv(i)) THEN
810	! cc nrlmd Introduction du taux de mortalité des poches et
811	! test sur sigmaw_max=0.4
812	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
813	IF (sigmaw(i)>=sigmaw_max) THEN
814	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
815	ELSE
816	death_rate(i) = 0.
817	END IF
818
819	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
820	dtimesub
821	! $ - nat_rate(i)sigmaw(i)dtimesub
822	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
823	! c $ death_rate(i),ktop(i),kupper(i)',
824	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
825	! c $ death_rate(i),ktop(i),kupper(i)
826
827	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
828	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
829	! wdens = wdens0/(10.*sigmaw)
830	! sigmaw =max(sigmaw,sigd_con)
831	! sigmaw =max(sigmaw,sigmad)
832	END IF
833	END DO
834
835	! calcul de la difference de vitesse verticale poche - zone non perturbee
836	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
837	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on
838	! definit
839	! IM 060208 au niveau k=1..?
840	DO k = 1, klev
841	DO i = 1, klon
842	IF (wk_adv(i)) THEN !!! nrlmd
843	dp_deltomg(i, k) = 0.
844	END IF
845	END DO
846	END DO
847	DO k = 1, klev + 1
848	DO i = 1, klon
849	IF (wk_adv(i)) THEN !!! nrlmd
850	omg(i, k) = 0.
851	END IF
852	END DO
853	END DO
854
855	DO i = 1, klon
856	IF (wk_adv(i)) THEN
857	z(i) = 0.
858	omg(i, 1) = 0.
859	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
860	END IF
861	END DO
862
863	DO k = 2, klev
864	DO i = 1, klon
865	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
866	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
867	z(i) = z(i) + dz(i)
868	dp_deltomg(i, k) = dp_deltomg(i, 1)
869	omg(i, k) = dp_deltomg(i, 1)*z(i)
870	END IF
871	END DO
872	END DO
873
874	DO i = 1, klon
875	IF (wk_adv(i)) THEN
876	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
877	ztop(i) = z(i) + dztop(i)
878	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
879	END IF
880	END DO
881
882	! -----------------
883	! From m/s to Pa/s
884	! -----------------
885
886	DO i = 1, klon
887	IF (wk_adv(i)) THEN
888	omgtop(i) = -rho(i, ktop(i))rgomgtop(i)
889	dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
890	END IF
891	END DO
892
893	DO k = 1, klev
894	DO i = 1, klon
895	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
896	omg(i, k) = -rho(i, k)rgomg(i, k)
897	dp_deltomg(i, k) = dp_deltomg(i, 1)
898	END IF
899	END DO
900	END DO
901
902	! raccordement lineaire de omg de ptop a pupper
903
904	DO i = 1, klon
905	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
906	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
907	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
908	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
909	(ptop(i)-pupper(i))
910	END IF
911	END DO
912
913	! c DO i=1,klon
914	! c print*,'Pente entre 0 et kupper (référence)'
915	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
916	! c print*,'Pente entre ktop et kupper'
917	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
918	! c ENDDO
919	! c
920	DO k = 1, klev
921	DO i = 1, klon
922	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
923	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
924	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
925	END IF
926	END DO
927	END DO
928	! cc nrlmd
929	! c DO i=1,klon
930	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
931	! c END DO
932	! cc
933
934
935	! -- Compute wake average vertical velocity omgbw
936
937
938	DO k = 1, klev + 1
939	DO i = 1, klon
940	IF (wk_adv(i)) THEN
941	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
942	END IF
943	END DO
944	END DO
945	! -- and its vertical gradient dp_omgbw
946
947	DO k = 1, klev
948	DO i = 1, klon
949	IF (wk_adv(i)) THEN
950	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
951	END IF
952	END DO
953	END DO
954
955	! -- Upstream coefficients for omgb velocity
956	! -- (alpha_up(k) is the coefficient of the value at level k)
957	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
958	DO k = 1, klev
959	DO i = 1, klon
960	IF (wk_adv(i)) THEN
961	alpha_up(i, k) = 0.
962	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
963	END IF
964	END DO
965	END DO
966
967	! Matrix expressing [The,deltatw] from [Th1,Th2]
968
969	DO i = 1, klon
970	IF (wk_adv(i)) THEN
971	rre1(i) = 1. - sigmaw(i)
972	rre2(i) = sigmaw(i)
973	END IF
974	END DO
975	rrd1 = -1.
976	rrd2 = 1.
977
978	! -- Get [Th1,Th2], dth and [q1,q2]
979
980	DO k = 1, klev
981	DO i = 1, klon
982	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
983	dth(i, k) = deltatw(i, k)/ppi(i, k)
984	th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
985	th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
986	q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
987	q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
988	END IF
989	END DO
990	END DO
991
992	DO i = 1, klon
993	IF (wk_adv(i)) THEN !!! nrlmd
994	d_th1(i, 1) = 0.
995	d_th2(i, 1) = 0.
996	d_dth(i, 1) = 0.
997	d_q1(i, 1) = 0.
998	d_q2(i, 1) = 0.
999	d_dq(i, 1) = 0.
1000	END IF
1001	END DO
1002
1003	DO k = 2, klev
1004	DO i = 1, klon
1005	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1006	d_th1(i, k) = th1(i, k-1) - th1(i, k)
1007	d_th2(i, k) = th2(i, k-1) - th2(i, k)
1008	d_dth(i, k) = dth(i, k-1) - dth(i, k)
1009	d_q1(i, k) = q1(i, k-1) - q1(i, k)
1010	d_q2(i, k) = q2(i, k-1) - q2(i, k)
1011	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
1012	END IF
1013	END DO
1014	END DO
1015
1016	DO i = 1, klon
1017	IF (wk_adv(i)) THEN
1018	omgbdth(i, 1) = 0.
1019	omgbdq(i, 1) = 0.
1020	END IF
1021	END DO
1022
1023	DO k = 2, klev
1024	DO i = 1, klon
1025	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
1026	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
1027	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
1028	END IF
1029	END DO
1030	END DO
1031
1032	! -----------------------------------------------------------------
1033	DO k = 1, klev
1034	DO i = 1, klon
1035	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1036	! -----------------------------------------------------------------
1037
1038	! Compute redistribution (advective) term
1039
1040	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1041	(rrd1omg(i,k)sigmaw(i)*d_th1(i,k) - &
1042	rrd2omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1)- &
1043	(1.-alpha_up(i,k))*omgbdth(i,k)- &
1044	alpha_up(i,k+1)omgbdth(i,k+1))ppi(i, k)
1045	! print*,'d_deltatw=',d_deltatw(i,k)
1046
1047	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1048	(rrd1omg(i,k)sigmaw(i)*d_q1(i,k)- &
1049	rrd2omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1)- &
1050	(1.-alpha_up(i,k))*omgbdq(i,k)- &
1051	alpha_up(i,k+1)*omgbdq(i,k+1))
1052	! print*,'d_deltaqw=',d_deltaqw(i,k)
1053
1054	! and increment large scale tendencies
1055
1056
1057
1058
1059	! C
1060	! -----------------------------------------------------------------
1061	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)- &
1062	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/ &
1063	(ph(i,k)-ph(i,k+1)) &
1064	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
1065	(ph(i,k)-ph(i,k+1)) )*ppi(i, k)
1066
1067	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)- &
1068	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/ &
1069	(ph(i,k)-ph(i,k+1)) &
1070	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
1071	(ph(i,k)-ph(i,k+1)) )
1072	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1073	d_te(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
1074
1075	d_qe(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
1076
1077	END IF
1078	! cc
1079	END DO
1080	END DO
1081	! ------------------------------------------------------------------
1082
1083	! Increment state variables
1084
1085	DO k = 1, klev
1086	DO i = 1, klon
1087	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1088	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1089	! cc
1090
1091
1092
1093	! Coefficient de répartition
1094
1095	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1096	(ph(i,kupper(i))-ph(i,1))
1097	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
1098	(p(i,1)-ph(i,kupper(i)))
1099
1100
1101	! Reintroduce compensating subsidence term.
1102
1103	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1104	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1105	! . /(1-sigmaw)
1106	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1107	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1108	! . /(1-sigmaw)
1109
1110	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
1111	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
1112	! . /(1-sigmaw)
1113	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
1114	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
1115	! . /(1-sigmaw)
1116
1117	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1118	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1119	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1120
1121	!
1122
1123	! cc nrlmd Prise en compte du taux de mortalité
1124	! cc Définitions de entr, detr
1125	detr(i, k) = 0.
1126
1127	entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1128	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1129
1130	spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1131	! cc spread(i,k) =
1132	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
1133	! cc $ sigmaw(i)
1134
1135
1136	! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
1137	! Jingmei
1138
1139	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1140	! & Tgw(i,k),deltatw(i,k)
1141	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
1142	dtimesub
1143	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1144	ff(i) = d_deltatw(i, k)/dtimesub
1145
1146	! Sans GW
1147
1148	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
1149
1150	! GW formule 1
1151
1152	! deltatw(k) = deltatw(k)+dtimesub*
1153	! $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
1154
1155	! GW formule 2
1156
1157	IF (dtimesub*tgw(i,k)<1.E-10) THEN
1158	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
1159	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1160	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
1161	tgw(i,k)*deltatw(i,k) )
1162	ELSE
1163	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i,k)))* &
1164	(ff(i)+dtke(i,k) - &
1165	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1166	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - &
1167	tgw(i,k)*deltatw(i,k) )
1168	END IF
1169
1170	dth(i, k) = deltatw(i, k)/ppi(i, k)
1171
1172	gg(i) = d_deltaqw(i, k)/dtimesub
1173
1174	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
1175	entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
1176	(death_rate(i)sigmaw(i)+detr(i,k))deltaqw(i,k)/(1.-sigmaw(i)))
1177	! cc
1178
1179	! cc nrlmd
1180	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1181	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1182	! cc
1183	END IF
1184	END DO
1185	END DO
1186
1187
1188	! Scale tendencies so that water vapour remains positive in w and x.
1189
1190	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
1191	d_deltaqw, sigmaw, d_sigmaw, alpha)
1192
1193	! cc nrlmd
1194	! c print*,'alpha'
1195	! c do i=1,klon
1196	! c print*,alpha(i)
1197	! c end do
1198	! cc
1199	DO k = 1, klev
1200	DO i = 1, klon
1201	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1202	d_te(i, k) = alpha(i)*d_te(i, k)
1203	d_qe(i, k) = alpha(i)*d_qe(i, k)
1204	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1205	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1206	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1207	END IF
1208	END DO
1209	END DO
1210	DO i = 1, klon
1211	IF (wk_adv(i)) THEN
1212	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1213	END IF
1214	END DO
1215
1216	! Update large scale variables and wake variables
1217	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1218	! IM 060208 DO k = 1,kupper(i)
1219	DO k = 1, klev
1220	DO i = 1, klon
1221	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1222	dtls(i, k) = dtls(i, k) + d_te(i, k)
1223	dqls(i, k) = dqls(i, k) + d_qe(i, k)
1224	! cc nrlmd
1225	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1226	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1227	! cc
1228	END IF
1229	END DO
1230	END DO
1231	DO k = 1, klev
1232	DO i = 1, klon
1233	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1234	te(i, k) = te0(i, k) + dtls(i, k)
1235	qe(i, k) = qe0(i, k) + dqls(i, k)
1236	the(i, k) = te(i, k)/ppi(i, k)
1237	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1238	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1239	dth(i, k) = deltatw(i, k)/ppi(i, k)
1240	! c print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
1241	! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1242	END IF
1243	END DO
1244	END DO
1245	DO i = 1, klon
1246	IF (wk_adv(i)) THEN
1247	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1248	!jyg<
1249	d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
1250	!>jyg
1251	END IF
1252	END DO
1253
1254
1255	! Determine Ptop from buoyancy integral
1256	! ---------------------------------------
1257
1258	! - 1/ Pressure of the level where dth changes sign.
1259
1260	DO i = 1, klon
1261	IF (wk_adv(i)) THEN
1262	ptop_provis(i) = ph(i, 1)
1263	END IF
1264	END DO
1265
1266	DO k = 2, klev
1267	DO i = 1, klon
1268	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1269	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1270	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1271	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1272	END IF
1273	END DO
1274	END DO
1275
1276	! - 2/ dth integral
1277
1278	DO i = 1, klon
1279	IF (wk_adv(i)) THEN !!! nrlmd
1280	sum_dth(i) = 0.
1281	dthmin(i) = -delta_t_min
1282	z(i) = 0.
1283	END IF
1284	END DO
1285
1286	DO k = 1, klev
1287	DO i = 1, klon
1288	IF (wk_adv(i)) THEN
1289	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1290	IF (dz(i)>0) THEN
1291	z(i) = z(i) + dz(i)
1292	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1293	dthmin(i) = amin1(dthmin(i), dth(i,k))
1294	END IF
1295	END IF
1296	END DO
1297	END DO
1298
1299	! - 3/ height of triangle with area= sum_dth and base = dthmin
1300
1301	DO i = 1, klon
1302	IF (wk_adv(i)) THEN
1303	hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1304	hw(i) = amax1(hwmin, hw(i))
1305	END IF
1306	END DO
1307
1308	! - 4/ now, get Ptop
1309
1310	DO i = 1, klon
1311	IF (wk_adv(i)) THEN !!! nrlmd
1312	ktop(i) = 0
1313	z(i) = 0.
1314	END IF
1315	END DO
1316
1317	DO k = 1, klev
1318	DO i = 1, klon
1319	IF (wk_adv(i)) THEN
1320	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
1321	IF (dz(i)>0) THEN
1322	z(i) = z(i) + dz(i)
1323	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
1324	ktop(i) = k
1325	END IF
1326	END IF
1327	END DO
1328	END DO
1329
1330	! 4.5/Correct ktop and ptop
1331
1332	DO i = 1, klon
1333	IF (wk_adv(i)) THEN
1334	ptop_new(i) = ptop(i)
1335	END IF
1336	END DO
1337
1338	DO k = klev, 2, -1
1339	DO i = 1, klon
1340	! IM v3JYG; IF (k .GE. ktop(i)
1341	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1342	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1343	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1344	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1345	END IF
1346	END DO
1347	END DO
1348
1349
1350	DO i = 1, klon
1351	IF (wk_adv(i)) THEN
1352	ptop(i) = ptop_new(i)
1353	END IF
1354	END DO
1355
1356	DO k = klev, 1, -1
1357	DO i = 1, klon
1358	IF (wk_adv(i)) THEN !!! nrlmd
1359	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
1360	END IF
1361	END DO
1362	END DO
1363
1364	! 5/ Set deltatw & deltaqw to 0 above kupper
1365
1366	DO k = 1, klev
1367	DO i = 1, klon
1368	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1369	deltatw(i, k) = 0.
1370	deltaqw(i, k) = 0.
1371	d_deltatw2(i,k) = -deltatw0(i,k)
1372	d_deltaqw2(i,k) = -deltaqw0(i,k)
1373	END IF
1374	END DO
1375	END DO
1376
1377
1378	! -------------Cstar computation---------------------------------
1379	DO i = 1, klon
1380	IF (wk_adv(i)) THEN !!! nrlmd
1381	sum_thu(i) = 0.
1382	sum_tu(i) = 0.
1383	sum_qu(i) = 0.
1384	sum_thvu(i) = 0.
1385	sum_dth(i) = 0.
1386	sum_dq(i) = 0.
1387	sum_rho(i) = 0.
1388	sum_dtdwn(i) = 0.
1389	sum_dqdwn(i) = 0.
1390
1391	av_thu(i) = 0.
1392	av_tu(i) = 0.
1393	av_qu(i) = 0.
1394	av_thvu(i) = 0.
1395	av_dth(i) = 0.
1396	av_dq(i) = 0.
1397	av_rho(i) = 0.
1398	av_dtdwn(i) = 0.
1399	av_dqdwn(i) = 0.
1400	END IF
1401	END DO
1402
1403	! Integrals (and wake top level number)
1404	! --------------------------------------
1405
1406	! Initialize sum_thvu to 1st level virt. pot. temp.
1407
1408	DO i = 1, klon
1409	IF (wk_adv(i)) THEN !!! nrlmd
1410	z(i) = 1.
1411	dz(i) = 1.
1412	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
1413	sum_dth(i) = 0.
1414	END IF
1415	END DO
1416
1417	DO k = 1, klev
1418	DO i = 1, klon
1419	IF (wk_adv(i)) THEN !!! nrlmd
1420	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1421	IF (dz(i)>0) THEN
1422	z(i) = z(i) + dz(i)
1423	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1424	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1425	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1426	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
1427	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1428	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1429	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1430	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1431	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1432	END IF
1433	END IF
1434	END DO
1435	END DO
1436
1437	DO i = 1, klon
1438	IF (wk_adv(i)) THEN !!! nrlmd
1439	hw0(i) = z(i)
1440	END IF
1441	END DO
1442
1443
1444	! - WAPE and mean forcing computation
1445	! ---------------------------------------
1446
1447	! ---------------------------------------
1448
1449	! Means
1450
1451	DO i = 1, klon
1452	IF (wk_adv(i)) THEN !!! nrlmd
1453	av_thu(i) = sum_thu(i)/hw0(i)
1454	av_tu(i) = sum_tu(i)/hw0(i)
1455	av_qu(i) = sum_qu(i)/hw0(i)
1456	av_thvu(i) = sum_thvu(i)/hw0(i)
1457	av_dth(i) = sum_dth(i)/hw0(i)
1458	av_dq(i) = sum_dq(i)/hw0(i)
1459	av_rho(i) = sum_rho(i)/hw0(i)
1460	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1461	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1462
1463	wape(i) = -rghw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
1464	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
1465	END IF
1466	END DO
1467
1468	! Filter out bad wakes
1469
1470	DO k = 1, klev
1471	DO i = 1, klon
1472	IF (wk_adv(i)) THEN !!! nrlmd
1473	IF (wape(i)<0.) THEN
1474	deltatw(i, k) = 0.
1475	deltaqw(i, k) = 0.
1476	dth(i, k) = 0.
1477	d_deltatw2(i,k) = -deltatw0(i,k)
1478	d_deltaqw2(i,k) = -deltaqw0(i,k)
1479	END IF
1480	END IF
1481	END DO
1482	END DO
1483
1484	DO i = 1, klon
1485	IF (wk_adv(i)) THEN !!! nrlmd
1486	IF (wape(i)<0.) THEN
1487	wape(i) = 0.
1488	cstar(i) = 0.
1489	hw(i) = hwmin
1490	!jyg<
1491	!! sigmaw(i) = max(sigmad, sigd_con(i))
1492	sigmaw_targ = max(sigmad, sigd_con(i))
1493	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1494	sigmaw(i) = sigmaw_targ
1495	!>jyg
1496	fip(i) = 0.
1497	gwake(i) = .FALSE.
1498	ELSE
1499	cstar(i) = starksqrt(2.wape(i))
1500	gwake(i) = .TRUE.
1501	END IF
1502	END IF
1503	END DO
1504
1505	END DO ! end sub-timestep loop
1506
1507
1508
1509	! ----------------------------------------------------------
1510	! Determine wake final state; recompute wape, cstar, ktop;
1511	! filter out bad wakes.
1512	! ----------------------------------------------------------
1513
1514	! 2.1 - Undisturbed area and Wake integrals
1515	! ---------------------------------------------------------
1516
1517	DO i = 1, klon
1518	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
1519	IF (ok_qx_qw(i)) THEN
1520	! cc
1521	z(i) = 0.
1522	sum_thu(i) = 0.
1523	sum_tu(i) = 0.
1524	sum_qu(i) = 0.
1525	sum_thvu(i) = 0.
1526	sum_dth(i) = 0.
1527	sum_dq(i) = 0.
1528	sum_rho(i) = 0.
1529	sum_dtdwn(i) = 0.
1530	sum_dqdwn(i) = 0.
1531
1532	av_thu(i) = 0.
1533	av_tu(i) = 0.
1534	av_qu(i) = 0.
1535	av_thvu(i) = 0.
1536	av_dth(i) = 0.
1537	av_dq(i) = 0.
1538	av_rho(i) = 0.
1539	av_dtdwn(i) = 0.
1540	av_dqdwn(i) = 0.
1541	END IF
1542	END DO
1543	! Potential temperatures and humidity
1544	! ----------------------------------------------------------
1545
1546	DO k = 1, klev
1547	DO i = 1, klon
1548	! cc nrlmd IF ( wk_adv(i)) THEN
1549	IF (ok_qx_qw(i)) THEN
1550	! cc
1551	rho(i, k) = p(i, k)/(rd*te(i,k))
1552	IF (k==1) THEN
1553	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
1554	zhh(i, k) = 0
1555	ELSE
1556	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
1557	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
1558	END IF
1559	the(i, k) = te(i, k)/ppi(i, k)
1560	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1561	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
1562	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
1563	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
1564	dth(i, k) = deltatw(i, k)/ppi(i, k)
1565	END IF
1566	END DO
1567	END DO
1568
1569	! Integrals (and wake top level number)
1570	! -----------------------------------------------------------
1571
1572	! Initialize sum_thvu to 1st level virt. pot. temp.
1573
1574	DO i = 1, klon
1575	! cc nrlmd IF ( wk_adv(i)) THEN
1576	IF (ok_qx_qw(i)) THEN
1577	! cc
1578	z(i) = 1.
1579	dz(i) = 1.
1580	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
1581	sum_dth(i) = 0.
1582	END IF
1583	END DO
1584
1585	DO k = 1, klev
1586	DO i = 1, klon
1587	! cc nrlmd IF ( wk_adv(i)) THEN
1588	IF (ok_qx_qw(i)) THEN
1589	! cc
1590	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1591	IF (dz(i)>0) THEN
1592	z(i) = z(i) + dz(i)
1593	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
1594	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
1595	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
1596	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
1597	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1598	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1599	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
1600	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1601	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1602	END IF
1603	END IF
1604	END DO
1605	END DO
1606
1607	DO i = 1, klon
1608	! cc nrlmd IF ( wk_adv(i)) THEN
1609	IF (ok_qx_qw(i)) THEN
1610	! cc
1611	hw0(i) = z(i)
1612	END IF
1613	END DO
1614
1615	! - WAPE and mean forcing computation
1616	! -------------------------------------------------------------
1617
1618	! Means
1619
1620	DO i = 1, klon
1621	! cc nrlmd IF ( wk_adv(i)) THEN
1622	IF (ok_qx_qw(i)) THEN
1623	! cc
1624	av_thu(i) = sum_thu(i)/hw0(i)
1625	av_tu(i) = sum_tu(i)/hw0(i)
1626	av_qu(i) = sum_qu(i)/hw0(i)
1627	av_thvu(i) = sum_thvu(i)/hw0(i)
1628	av_dth(i) = sum_dth(i)/hw0(i)
1629	av_dq(i) = sum_dq(i)/hw0(i)
1630	av_rho(i) = sum_rho(i)/hw0(i)
1631	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1632	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1633
1634	wape2(i) = -rghw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
1635	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
1636	END IF
1637	END DO
1638
1639
1640
1641	! Prognostic variable update
1642	! ------------------------------------------------------------
1643
1644	! Filter out bad wakes
1645
1646	DO k = 1, klev
1647	DO i = 1, klon
1648	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
1649	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
1650	! cc
1651	deltatw(i, k) = 0.
1652	deltaqw(i, k) = 0.
1653	dth(i, k) = 0.
1654	d_deltatw2(i,k) = -deltatw0(i,k)
1655	d_deltaqw2(i,k) = -deltaqw0(i,k)
1656	END IF
1657	END DO
1658	END DO
1659
1660
1661	DO i = 1, klon
1662	! cc nrlmd IF ( wk_adv(i)) THEN
1663	IF (ok_qx_qw(i)) THEN
1664	! cc
1665	IF (wape2(i)<0.) THEN
1666	wape2(i) = 0.
1667	cstar2(i) = 0.
1668	hw(i) = hwmin
1669	!jyg<
1670	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
1671	sigmaw_targ = max(sigmad, sigd_con(i))
1672	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1673	sigmaw(i) = sigmaw_targ
1674	!>jyg
1675	fip(i) = 0.
1676	gwake(i) = .FALSE.
1677	ELSE
1678	IF (prt_level>=10) PRINT *, 'wape2>0'
1679	cstar2(i) = starksqrt(2.wape2(i))
1680	gwake(i) = .TRUE.
1681	END IF
1682	END IF
1683	END DO
1684
1685	DO i = 1, klon
1686	! cc nrlmd IF ( wk_adv(i)) THEN
1687	IF (ok_qx_qw(i)) THEN
1688	! cc
1689	ktopw(i) = ktop(i)
1690	END IF
1691	END DO
1692
1693	DO i = 1, klon
1694	! cc nrlmd IF ( wk_adv(i)) THEN
1695	IF (ok_qx_qw(i)) THEN
1696	! cc
1697	IF (ktopw(i)>0 .AND. gwake(i)) THEN
1698
1699	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
1700	! cc heff = 600.
1701	! Utilisation de la hauteur hw
1702	! c heff = 0.7*hw
1703	heff(i) = hw(i)
1704
1705	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
1706	sqrt(sigmaw(i)wdens(i)3.14)
1707	fip(i) = alpk*fip(i)
1708	! jyg2
1709	ELSE
1710	fip(i) = 0.
1711	END IF
1712	END IF
1713	END DO
1714
1715	! Limitation de sigmaw
1716
1717	! cc nrlmd
1718	! DO i=1,klon
1719	! IF (OK_qx_qw(i)) THEN
1720	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
1721	! ENDIF
1722	! ENDDO
1723	! cc
1724	DO k = 1, klev
1725	DO i = 1, klon
1726
1727	! cc nrlmd On maintient désormais constant sigmaw en régime
1728	! permanent
1729	! cc IF ((sigmaw(i).GT.sigmaw_max).or.
1730	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1731	.NOT. ok_qx_qw(i)) THEN
1732	! cc
1733	dtls(i, k) = 0.
1734	dqls(i, k) = 0.
1735	deltatw(i, k) = 0.
1736	deltaqw(i, k) = 0.
1737	d_deltatw2(i,k) = -deltatw0(i,k)
1738	d_deltaqw2(i,k) = -deltaqw0(i,k)
1739	END IF
1740	END DO
1741	END DO
1742
1743	! cc nrlmd On maintient désormais constant sigmaw en régime permanent
1744	DO i = 1, klon
1745	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
1746	.NOT. ok_qx_qw(i)) THEN
1747	ktopw(i) = 0
1748	wape(i) = 0.
1749	cstar(i) = 0.
1750	!!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
1751	!! hw(i) = hwmin !jyg
1752	!! sigmaw(i) = sigmad !jyg
1753	hw(i) = 0. !jyg
1754	sigmaw(i) = 0. !jyg
1755	fip(i) = 0.
1756	ELSE
1757	wape(i) = wape2(i)
1758	cstar(i) = cstar2(i)
1759	END IF
1760	! c print*,'wape wape2 ktopw OK_qx_qw =',
1761	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
1762	END DO
1763
1764	! -----------------------------------------------------------------
1765	! Get back to tendencies per second
1766
1767	DO k = 1, klev
1768	DO i = 1, klon
1769
1770	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1771	IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
1772	! cc
1773	dtls(i, k) = dtls(i, k)/dtime
1774	dqls(i, k) = dqls(i, k)/dtime
1775	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
1776	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
1777	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
1778	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
1779	! c $ ,death_rate(i)*sigmaw(i)
1780	END IF
1781	END DO
1782	END DO
1783	!jyg<
1784	DO i = 1, klon
1785	d_sigmaw2(i) = d_sigmaw2(i)/dtime
1786	d_wdens2(i) = d_wdens2(i)/dtime
1787	ENDDO
1788	!>jyg
1789
1790
1791
1792	RETURN
1793	END SUBROUTINE wake
1794
1795	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
1796	d_deltaqw, sigmaw, d_sigmaw, alpha)
1797	! ------------------------------------------------------
1798	! Dtermination du coefficient alpha tel que les tendances
1799	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
1800	! a une humidite positive dans la zone (x) et dans la zone (w).
1801	! ------------------------------------------------------
1802	IMPLICIT NONE
1803
1804	! Input
1805	REAL qe(nlon, nl), d_qe(nlon, nl)
1806	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
1807	REAL sigmaw(nlon), d_sigmaw(nlon)
1808	LOGICAL wk_adv(nlon)
1809	INTEGER nl, nlon
1810	! Output
1811	REAL alpha(nlon)
1812	! Internal variables
1813	REAL zeta(nlon, nl)
1814	REAL alpha1(nlon)
1815	REAL x, a, b, c, discrim
1816	REAL epsilon
1817	! DATA epsilon/1.e-15/
1818	INTEGER i,k
1819
1820	DO k = 1, nl
1821	DO i = 1, nlon
1822	IF (wk_adv(i)) THEN
1823	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
1824	zeta(i, k) = 0.
1825	ELSE
1826	zeta(i, k) = 1.
1827	END IF
1828	END IF
1829	END DO
1830	DO i = 1, nlon
1831	IF (wk_adv(i)) THEN
1832	x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
1833	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i) &
1834	(deltaqw(i,k)+d_deltaqw(i,k))
1835	a = -d_sigmaw(i)*d_deltaqw(i, k)
1836	b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
1837	deltaqw(i, k)*d_sigmaw(i)
1838	c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
1839	discrim = bb - 4.a*c
1840	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
1841	IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
1842	alpha1(i) = 1.
1843	ELSE
1844	IF (x>=0.) THEN
1845	alpha1(i) = 1.
1846	ELSE
1847	IF (a>0.) THEN
1848	alpha1(i) = 0.9min( (2.c)/(-b+sqrt(discrim)), &
1849	(-b+sqrt(discrim))/(2.*a) )
1850	ELSE IF (a==0.) THEN
1851	alpha1(i) = 0.9*(-c/b)
1852	ELSE
1853	! print*,'a,b,c discrim',a,b,c discrim
1854	alpha1(i) = 0.9max( (2.c)/(-b+sqrt(discrim)), &
1855	(-b+sqrt(discrim))/(2.*a))
1856	END IF
1857	END IF
1858	END IF
1859	alpha(i) = min(alpha(i), alpha1(i))
1860	END IF
1861	END DO
1862	END DO
1863
1864	RETURN
1865	END SUBROUTINE wake_vec_modulation
1866
1867
1868
1869

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

Download in other formats:

Original Format