Context Navigation

wake.F90 @ 5225

Last change on this file since 5225 was 3899, checked in by lguez, 4 years ago
Old work on cold pool death.
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: 64.6 KB

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

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

Download in other formats:

Original Format