Context Navigation

wake.F90 @ 5225

Last change on this file since 5225 was 2925, checked in by fcheruy, 7 years ago
Update tree branch to trunk version
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: 62.4 KB

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

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

Download in other formats:

Original Format