Context Navigation

lmdz_wake.f90 @ 5503

Last change on this file since 5503 was 5499, checked in by yann meurdesoif, 2 days ago
Adapt calwake and wake for automatic GPU port. YM
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: 120.4 KB

Line
1	MODULE lmdz_wake
2
3	USE lmdz_wake_ini, ONLY: CPPKEY_IOPHYS_WK
4
5	IMPLICIT NONE; PRIVATE
6
7	LOGICAL, PARAMETER :: phys_sub=.false.
8	LOGICAL :: first_call=.true.
9	!$OMP THREADPRIVATE(first_call)
10
11	PUBLIC wake, wake_first
12
13	CONTAINS
14
15	SUBROUTINE wake_first(klev, dtime)
16	USE lmdz_wake_ini , ONLY : wk_nsub
17	IMPLICIT NONE
18	INTEGER, INTENT(IN) :: klev
19	REAL, INTENT(IN) :: dtime
20	REAL :: dtimesub
21
22	dtimesub = dtime/wk_nsub
23	!
24	IF (first_call) THEN
25	IF (CPPKEY_IOPHYS_WK) THEN
26	IF (phys_sub) THEN
27	call iophys_ini(dtimesub,klev)
28	ELSE
29	call iophys_ini(dtime,klev)
30	ENDIF
31	END IF
32	first_call = .false.
33	ENDIF !(first_call)
34
35	END SUBROUTINE wake_first
36
37	SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, &
38	tb0, qb0, omgb, &
39	dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
40	sigd_con, Cin, &
41	deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, & ! state variables
42	dth, hw, wape, fip, gfl, &
43	dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
44	dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
45	d_deltat_gw, & ! tendencies
46	d_deltatw2, d_deltaqw2, d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2) ! tendencies
47
48
49	! **************************************************************
50	! *
51	! WAKE *
52	! retour a un Pupper fixe *
53	! *
54	! written by : GRANDPEIX Jean-Yves 09/03/2000 *
55	! modified by : ROEHRIG Romain 01/29/2007 *
56	! **************************************************************
57
58
59	USE lmdz_wake_ini , ONLY : wake_ini
60	USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD
61	USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper
62	USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
63	USE lmdz_wake_ini , ONLY : ok_bug_gfl
64	USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensinit, wdensthreshold
65	USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
66	USE lmdz_wake_ini , ONLY : iflag_wk_profile
67	USE lmdz_wake_ini , ONLY : smallestreal,wk_nsub
68
69
70	IMPLICIT NONE
71	! ============================================================================
72
73
74	! But : Decrire le comportement des poches froides apparaissant dans les
75	! grands systemes convectifs, et fournir l'energie disponible pour
76	! le declenchement de nouvelles colonnes convectives.
77
78	! State variables :
79	! deltatw : temperature difference between wake and off-wake regions
80	! deltaqw : specific humidity difference between wake and off-wake regions
81	! sigmaw : fractional area covered by wakes.
82	! asigmaw : fractional area covered by active wakes.
83	! wdens : number of wakes per unit area
84	! awdens : number of active wakes per unit area
85
86	! Variable de sortie :
87
88	! wape : WAke Potential Energy
89	! fip : Front Incident Power (W/m2) - ALP
90	! gfl : Gust Front Length per unit area (m-1)
91	! dtls : large scale temperature tendency due to wake
92	! dqls : large scale humidity tendency due to wake
93	! hw : wake top hight (given by hw*deltatw(1)/2=wape)
94	! dp_omgb : vertical gradient of large scale omega
95	! awdens : densite de poches actives
96	! wdens : densite de poches
97	! omgbdth: flux of Delta_Theta transported by LS omega
98	! dtKE : differential heating (wake - unpertubed)
99	! dqKE : differential moistening (wake - unpertubed)
100	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
101	! dp_deltomg : vertical gradient of omg (s-1)
102	! wkspread : spreading term in d_t_wake and d_q_wake
103	! deltatw : updated temperature difference (T_w-T_u).
104	! deltaqw : updated humidity difference (q_w-q_u).
105	! sigmaw : updated wake fractional area.
106	! asigmaw : updated active wake fractional area.
107	! d_deltat_gw : delta T tendency due to GW
108
109	! Variables d'entree :
110
111	! aire : aire de la maille
112	! tb0 : horizontal average of temperature (K)
113	! qb0 : horizontal average of humidity (kg/kg)
114	! omgb : vitesse verticale moyenne sur la maille (Pa/s)
115	! dtdwn: source de chaleur due aux descentes (K/s)
116	! dqdwn: source d'humidite due aux descentes (kg/kg/s)
117	! dta : source de chaleur due courants satures et detrain (K/s)
118	! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
119	! wgen : number of wakes generated per unit area and per sec (/m^2/s)
120	! amdwn: flux de masse total des descentes, par unite de
121	! surface de la maille (kg/m2/s)
122	! amup : flux de masse total des ascendances, par unite de
123	! surface de la maille (kg/m2/s)
124	! sigd_con:
125	! Cin : convective inhibition
126	! p : pressions aux milieux des couches (Pa)
127	! ph : pressions aux interfaces (Pa)
128	! pi : (p/p_0)**kapa (adim)
129	! dtime: increment temporel (s)
130
131	! Variables internes :
132
133	! rho : mean density at P levels
134	! rhoh : mean density at Ph levels
135	! tb : mean temperature \| may change within
136	! qb : mean humidity \| sub-time-stepping
137	! thb : mean potential temperature
138	! thx : potential temperature in (x) area
139	! tx : temperature in (x) area
140	! qx : humidity in (x) area
141	! dp_omgb: vertical gradient og LS omega
142	! omgbw : wake average vertical omega
143	! dp_omgbw: vertical gradient of omgbw
144	! omgbdq : flux of Delta_q transported by LS omega
145	! dth : potential temperature diff. wake-undist.
146	! th1 : first pot. temp. for vertical advection (=thx)
147	! th2 : second pot. temp. for vertical advection (=thw)
148	! q1 : first humidity for vertical advection
149	! q2 : second humidity for vertical advection
150	! d_deltatw : redistribution term for deltatw
151	! d_deltaqw : redistribution term for deltaqw
152	! deltatw0 : initial deltatw
153	! deltaqw0 : initial deltaqw
154	! hw0 : wake top hight (defined as the altitude at which deltatw=0)
155	! amflux : horizontal mass flux through wake boundary
156	! wdens_ref: initial number of wakes per unit area (3D) or per
157	! unit length (2D), at the beginning of each time step
158	! Tgw : 1 sur la periode de onde de gravite
159	! Cgw : vitesse de propagation de onde de gravite
160	! LL : distance between 2 wakes
161
162	! -------------------------------------------------------------------------
163	! Declaration de variables
164	! -------------------------------------------------------------------------
165
166
167	! Arguments en entree
168	! --------------------
169
170	INTEGER, INTENT(IN) :: klon,klev
171	INTEGER, DIMENSION (klon), INTENT(IN) :: znatsurf
172	REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi
173	REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph
174	REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb
175	REAL, INTENT(IN) :: dtime
176	REAL, DIMENSION (klon, klev), INTENT(IN) :: tb0, qb0
177	REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn
178	REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup
179	REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa
180	REAL, DIMENSION (klon), INTENT(IN) :: wgen
181	REAL, DIMENSION (klon), INTENT(IN) :: sigd_con
182	REAL, DIMENSION (klon), INTENT(IN) :: Cin
183
184	!
185	! Input/Output
186	! State variables
187	REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw
188	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw
189	REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw
190	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens
191	REAL, DIMENSION (klon), INTENT(INOUT) :: awdens
192
193	! Sorties
194	! --------
195
196	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth
197	REAL, DIMENSION (klon, klev), INTENT(OUT) :: tx, qx
198	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls
199	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke
200	REAL, DIMENSION (klon, klev), INTENT(OUT) :: wkspread ! unused (jyg)
201	REAL, DIMENSION (klon, klev), INTENT(OUT) :: omgbdth, omg
202	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg
203	REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar
204	INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw
205	! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
206	! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
207	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw
208	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2
209	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw2, d_asigmaw2, d_wdens2, d_awdens2
210
211	! Variables internes
212	! -------------------
213
214	! Variables a fixer
215
216	REAL :: delta_t_min
217	REAL :: dtimesub
218	REAL :: wdens0
219	! IM 080208
220	LOGICAL, DIMENSION (klon) :: gwake
221
222	! Variables de sauvegarde
223	REAL, DIMENSION (klon, klev) :: deltatw0
224	REAL, DIMENSION (klon, klev) :: deltaqw0
225	REAL, DIMENSION (klon, klev) :: tb, qb
226
227	! Variables liees a la dynamique de population 1
228	REAL, DIMENSION(klon) :: act
229	REAL, DIMENSION(klon) :: rad_wk, tau_wk_inv
230	REAL, DIMENSION(klon) :: f_shear
231	REAL, DIMENSION(klon) :: drdt
232
233	! Variables liees a la dynamique de population 2
234	REAL, DIMENSION(klon) :: cont_fact
235
236	! Variables liees a la dynamique de population 3
237	REAL, DIMENSION(klon) :: arad_wk, irad_wk
238
239	!! REAL, DIMENSION(klon) :: d_sig_gen, d_sig_death, d_sig_col
240	REAL, DIMENSION(klon) :: wape1_act, wape2_act
241	LOGICAL, DIMENSION (klon) :: kill_wake
242	REAL :: drdt_pos
243	REAL :: tau_wk_inv_min
244	! Some components of the tendencies of state variables
245	REAL, DIMENSION (klon) :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_spread2, d_sig_bnd2
246	REAL, DIMENSION (klon) :: d_asig_death2, d_asig_aicol2, d_asig_iicol2, d_asig_spread2, d_asig_bnd2
247	REAL, DIMENSION (klon) :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
248	REAL, DIMENSION (klon) :: d_adens_death2, d_adens_icol2, d_adens_acol2, d_adens_bnd2
249
250	! Variables pour les GW
251	REAL, DIMENSION (klon) :: ll
252	REAL, DIMENSION (klon, klev) :: n2
253	REAL, DIMENSION (klon, klev) :: cgw
254	REAL, DIMENSION (klon, klev) :: tgw
255
256	! Variables liees au calcul de hw
257	REAL, DIMENSION (klon) :: ptop
258	REAL, DIMENSION (klon) :: sum_dth
259	REAL, DIMENSION (klon) :: dthmin
260	REAL, DIMENSION (klon) :: z, dz, hw0
261	INTEGER, DIMENSION (klon) :: ktop, kupper
262
263	! Variables liees au test de la forme triangulaire du profil de Delta_theta
264	REAL, DIMENSION (klon) :: sum_half_dth
265	REAL, DIMENSION (klon) :: dz_half
266
267	! Sub-timestep tendencies and related variables
268	REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw
269	REAL, DIMENSION (klon, klev) :: d_tb, d_qb
270	REAL, DIMENSION (klon) :: d_wdens, d_awdens, d_sigmaw, d_asigmaw
271	REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
272	REAL, DIMENSION (klon) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
273	REAL, DIMENSION (klon) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
274	REAL, DIMENSION (klon) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
275	REAL, DIMENSION (klon) :: agfl !! gust front length of active wakes
276	!! per unit area
277	REAL, DIMENSION (klon) :: alpha, alpha_tot
278	REAL, DIMENSION (klon) :: q0_min, q1_min
279	LOGICAL, DIMENSION (klon) :: wk_adv, ok_qx_qw
280
281	! Autres variables internes
282	INTEGER ::isubstep, k, i, igout
283
284	REAL :: wdensmin
285
286	REAL :: sigmaw_targ
287	REAL :: wdens_targ
288	REAL :: d_sigmaw_targ
289	REAL :: d_wdens_targ
290
291	REAL, DIMENSION (klon) :: sum_thx, sum_tx, sum_qx, sum_thvx
292	REAL, DIMENSION (klon) :: sum_dq
293	REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
294	REAL, DIMENSION (klon) :: av_thx, av_tx, av_qx, av_thvx
295	REAL, DIMENSION (klon) :: av_dth, av_dq
296	REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
297
298	REAL, DIMENSION (klon, klev) :: rho
299	REAL, DIMENSION (klon, klev+1) :: rhoh
300	REAL, DIMENSION (klon, klev) :: zh
301	REAL, DIMENSION (klon, klev+1) :: zhh
302
303	REAL, DIMENSION (klon, klev) :: thb, thx
304
305	REAL, DIMENSION (klon, klev) :: omgbw
306	REAL, DIMENSION (klon) :: pupper
307	REAL, DIMENSION (klon) :: omgtop
308	REAL, DIMENSION (klon, klev) :: dp_omgbw
309	REAL, DIMENSION (klon) :: ztop, dztop
310	REAL, DIMENSION (klon, klev) :: alpha_up
311
312	REAL, DIMENSION (klon) :: rre1, rre2
313	REAL :: rrd1, rrd2
314	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
315	REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
316	REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
317	REAL, DIMENSION (klon, klev) :: omgbdq
318
319	REAL, DIMENSION (klon) :: ff, gg
320	REAL, DIMENSION (klon) :: wape2, cstar2, heff
321
322	REAL, DIMENSION (klon, klev) :: crep
323
324	REAL, DIMENSION (klon, klev) :: ppi
325
326	! cc nrlmd
327	REAL, DIMENSION (klon) :: death_rate
328	!! REAL, DIMENSION (klon) :: nat_rate
329	REAL, DIMENSION (klon, klev) :: entr
330	REAL, DIMENSION (klon, klev) :: detr
331
332	REAL, DIMENSION(klon) :: sigmaw_in, asigmaw_in ! pour les prints
333	REAL, DIMENSION(klon) :: wdens_in, awdens_in ! pour les prints
334
335	!!-- variables liees au nouveau calcul de ptop et hw
336	REAL, DIMENSION (klon, klev) :: int_dth
337	REAL, DIMENSION (klon, klev) :: zzz, dzzz
338	REAL :: epsil
339	REAL, DIMENSION (klon) :: ptop1
340	INTEGER, DIMENSION (klon) :: ktop1
341	REAL, DIMENSION (klon) :: omega
342	REAL, DIMENSION (klon) :: h_zzz
343
344	!print*,'WAKE LJYFz'
345
346	! -------------------------------------------------------------------------
347	! Initialisations
348	! -------------------------------------------------------------------------
349	! ALON = 3.e5
350	! alon = 1.E6
351
352	! Provisionnal; to be suppressed when f_shear is parameterized
353	f_shear(:) = 1. ! 0. for strong shear, 1. for weak shear
354
355	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
356
357	! coefgw : Coefficient pour les ondes de gravite
358	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
359	! wdens : Densite surfacique de poche froide
360	! -------------------------------------------------------------------------
361
362	! cc nrlmd coefgw=10
363	! coefgw=1
364	! wdens0 = 1.0/(alon**2)
365	! cc nrlmd wdens = 1.0/(alon**2)
366	! cc nrlmd stark = 0.50
367	! CRtest
368	! cc nrlmd alpk=0.1
369	! alpk = 1.0
370	! alpk = 0.5
371	! alpk = 0.05
372	!
373	igout = klon/2+1/klon
374	!
375	! sub-time-stepping parameters
376	dtimesub = dtime/wk_nsub
377	!
378	IF (iflag_wk_pop_dyn == 0) THEN
379	! Initialisation de toutes des densites a wdens_ref.
380	! Les densites peuvent evoluer si les poches debordent
381	! (voir au tout debut de la boucle sur les substeps)
382	!jyg<
383	!! wdens(:) = wdens_ref
384	DO i = 1,klon
385	wdens(i) = wdens_ref(znatsurf(i)+1)
386	ENDDO
387	!>jyg
388	ENDIF ! (iflag_wk_pop_dyn == 0)
389	!
390	IF (iflag_wk_pop_dyn >=1) THEN
391	IF (iflag_wk_pop_dyn == 3) THEN
392	wdensmin = wdensthreshold
393	ELSE
394	wdensmin = wdensinit
395	ENDIF
396	ENDIF
397
398	! print*,'stark',stark
399	! print*,'alpk',alpk
400	! print*,'wdens',wdens
401	! print*,'coefgw',coefgw
402	! cc
403	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
404	! -------------------------------------------------------------------------
405
406	delta_t_min = 0.2
407
408	! 1. - Save initial values, initialize tendencies, initialize output fields
409	! ------------------------------------------------------------------------
410
411	!jyg<
412	!! DO k = 1, klev
413	!! DO i = 1, klon
414	!! ppi(i, k) = pi(i, k)
415	!! deltatw0(i, k) = deltatw(i, k)
416	!! deltaqw0(i, k) = deltaqw(i, k)
417	!! tb(i, k) = tb0(i, k)
418	!! qb(i, k) = qb0(i, k)
419	!! dtls(i, k) = 0.
420	!! dqls(i, k) = 0.
421	!! d_deltat_gw(i, k) = 0.
422	!! d_tb(i, k) = 0.
423	!! d_qb(i, k) = 0.
424	!! d_deltatw(i, k) = 0.
425	!! d_deltaqw(i, k) = 0.
426	!! ! IM 060508 beg
427	!! d_deltatw2(i, k) = 0.
428	!! d_deltaqw2(i, k) = 0.
429	!! ! IM 060508 end
430	!! END DO
431	!! END DO
432	ppi(:,:) = pi(:,:)
433	deltatw0(:,:) = deltatw(:,:)
434	deltaqw0(:,:) = deltaqw(:,:)
435	tb(:,:) = tb0(:,:)
436	qb(:,:) = qb0(:,:)
437	dtls(:,:) = 0.
438	dqls(:,:) = 0.
439	d_deltat_gw(:,:) = 0.
440	d_tb(:,:) = 0.
441	d_qb(:,:) = 0.
442	d_deltatw(:,:) = 0.
443	d_deltaqw(:,:) = 0.
444	d_deltatw2(:,:) = 0.
445	d_deltaqw2(:,:) = 0.
446
447	d_sig_gen2(:) = 0.
448	d_sig_death2(:) = 0.
449	d_sig_col2(:) = 0.
450	d_sig_spread2(:)= 0.
451	d_asig_death2(:) = 0.
452	d_asig_iicol2(:) = 0.
453	d_asig_aicol2(:) = 0.
454	d_asig_spread2(:)= 0.
455	d_asig_bnd2(:) = 0.
456	d_asigmaw2(:) = 0.
457	!
458	d_dens_gen2(:) = 0.
459	d_dens_death2(:) = 0.
460	d_dens_col2(:) = 0.
461	d_dens_bnd2(:) = 0.
462	d_wdens2(:) = 0.
463	d_adens_bnd2(:) = 0.
464	d_awdens2(:) = 0.
465	d_adens_death2(:) = 0.
466	d_adens_icol2(:) = 0.
467	d_adens_acol2(:) = 0.
468
469	IF (iflag_wk_act == 0) THEN
470	act(:) = 0.
471	ELSEIF (iflag_wk_act == 1) THEN
472	act(:) = 1.
473	ENDIF
474
475	!! DO i = 1, klon
476	!! sigmaw_in(i) = sigmaw(i)
477	!! END DO
478	sigmaw_in(:) = sigmaw(:)
479	asigmaw_in(:) = asigmaw(:)
480	!>jyg
481	!
482	IF (iflag_wk_pop_dyn >= 1) THEN
483	awdens_in(:) = awdens(:)
484	wdens_in(:) = wdens(:)
485	!! wdens(:) = wdens(:) + wgen(:)*dtime
486	!! d_wdens2(:) = wgen(:)*dtime
487	!! ELSE
488	ENDIF ! (iflag_wk_pop_dyn >= 1)
489
490
491	! sigmaw1=sigmaw
492	! IF (sigd_con.GT.sigmaw1) THEN
493	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
494	! ENDIF
495	IF (iflag_wk_pop_dyn >= 1) THEN
496	DO i = 1, klon
497	d_dens_gen2(i) = 0.
498	d_dens_death2(i) = 0.
499	d_dens_col2(i) = 0.
500	d_awdens2(i) = 0.
501	IF (wdens(i) < wdensthreshold) THEN
502	!! wdens_targ = max(wdens(i),wdensmin)
503	wdens_targ = max(wdens(i),wdensinit)
504	d_dens_bnd2(i) = wdens_targ - wdens(i)
505	d_wdens2(i) = wdens_targ - wdens(i)
506	wdens(i) = wdens_targ
507	ELSE
508	d_dens_bnd2(i) = 0.
509	d_wdens2(i) = 0.
510	ENDIF !! (wdens(i) < wdensthreshold)
511	END DO
512	IF (iflag_wk_pop_dyn >= 2) THEN
513	DO i = 1, klon
514	IF (awdens(i) < wdensthreshold) THEN
515	!! wdens_targ = min(max(awdens(i),wdensmin),wdens(i))
516	wdens_targ = min(max(awdens(i),wdensinit),wdens(i))
517	d_adens_bnd2(i) = wdens_targ - awdens(i)
518	d_awdens2(i) = wdens_targ - awdens(i)
519	awdens(i) = wdens_targ
520	ELSE
521	wdens_targ = min(awdens(i), wdens(i))
522	d_adens_bnd2(i) = wdens_targ - awdens(i)
523	d_awdens2(i) = wdens_targ - awdens(i)
524	awdens(i) = wdens_targ
525	ENDIF
526	END DO
527	ENDIF ! (iflag_wk_pop_dyn >= 2)
528	ELSE
529	DO i = 1, klon
530	d_awdens2(i) = 0.
531	d_wdens2(i) = 0.
532	END DO
533	ENDIF ! (iflag_wk_pop_dyn >= 1)
534	!
535	DO i = 1, klon
536	sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
537	d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
538	d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
539	sigmaw(i) = sigmaw_targ
540	END DO
541	!
542	IF (iflag_wk_pop_dyn == 3) THEN
543	DO i = 1, klon
544	IF ((wdens(i)-awdens(i)) <= smallestreal) THEN
545	sigmaw_targ = sigmaw(i)
546	ELSE
547	sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
548	ENDIF
549	d_asig_bnd2(i) = sigmaw_targ - asigmaw(i)
550	d_asigmaw2(i) = sigmaw_targ - asigmaw(i)
551	asigmaw(i) = sigmaw_targ
552	END DO
553	ENDIF ! (iflag_wk_pop_dyn == 3)
554
555	wape(:) = 0.
556	wape2(:) = 0.
557	d_sigmaw(:) = 0.
558	d_asigmaw(:) = 0.
559	ktopw(:) = 0
560	!
561	!<jyg
562	dth(:,:) = 0.
563	tx(:,:) = 0.
564	qx(:,:) = 0.
565	dtke(:,:) = 0.
566	dqke(:,:) = 0.
567	wkspread(:,:) = 0.
568	omgbdth(:,:) = 0.
569	omg(:,:) = 0.
570	dp_omgb(:,:) = 0.
571	dp_deltomg(:,:) = 0.
572	hw(:) = 0.
573	wape(:) = 0.
574	fip(:) = 0.
575	gfl(:) = 0.
576	cstar(:) = 0.
577	ktopw(:) = 0
578	!
579	! Vertical advection local variables
580	omgbw(:,:) = 0.
581	omgtop(:) = 0
582	dp_omgbw(:,:) = 0.
583	omgbdq(:,:) = 0.
584
585	!>jyg
586	!
587	IF (prt_level>=10) THEN
588	PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
589	PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
590	PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
591	PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
592	PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
593	PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
594	PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
595	PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
596	PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
597	ENDIF
598
599	! 2. - Prognostic part
600	! --------------------
601
602
603	! 2.1 - Undisturbed area and Wake integrals
604	! ---------------------------------------------------------
605
606	DO i = 1, klon
607	z(i) = 0.
608	ktop(i) = 0
609	kupper(i) = 0
610	sum_thx(i) = 0.
611	sum_tx(i) = 0.
612	sum_qx(i) = 0.
613	sum_thvx(i) = 0.
614	sum_dth(i) = 0.
615	sum_dq(i) = 0.
616	sum_dtdwn(i) = 0.
617	sum_dqdwn(i) = 0.
618
619	av_thx(i) = 0.
620	av_tx(i) = 0.
621	av_qx(i) = 0.
622	av_thvx(i) = 0.
623	av_dth(i) = 0.
624	av_dq(i) = 0.
625	av_dtdwn(i) = 0.
626	av_dqdwn(i) = 0.
627	END DO
628
629	! Distance between wakes
630	DO i = 1, klon
631	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
632	END DO
633	! Potential temperatures and humidity
634	! ----------------------------------------------------------
635	DO k = 1, klev
636	DO i = 1, klon
637	! write(,)'wake 1',i,k,RD,tb(i,k)
638	rho(i, k) = p(i, k)/(RD*tb(i,k))
639	! write(,)'wake 2',rho(i,k)
640	IF (k==1) THEN
641	! write(,)'wake 3',i,k,rd,tb(i,k)
642	rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
643	! write(,)'wake 4',i,k,rd,tb(i,k)
644	zhh(i, k) = 0
645	ELSE
646	! write(,)'wake 5',rd,(tb(i,k)+tb(i,k-1))
647	rhoh(i, k) = ph(i, k)2./(RD(tb(i,k)+tb(i,k-1)))
648	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
649	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
650	END IF
651	! write(,)'wake 7',ppi(i,k)
652	thb(i, k) = tb(i, k)/ppi(i, k)
653	thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
654	tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
655	qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
656	! write(,)'wake 8',(RD*(tb(i,k)+deltatw(i,k)))
657	dth(i, k) = deltatw(i, k)/ppi(i, k)
658	END DO
659	END DO
660
661	DO k = 1, klev - 1
662	DO i = 1, klon
663	IF (k==1) THEN
664	n2(i, k) = 0
665	ELSE
666	n2(i, k) = amax1(0., -RG*2/thb(i,k)rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ &
667	(p(i,k+1)-p(i,k-1)))
668	END IF
669	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
670
671	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
672	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
673	END DO
674	END DO
675
676	DO i = 1, klon
677	n2(i, klev) = 0
678	zh(i, klev) = 0
679	cgw(i, klev) = 0
680	tgw(i, klev) = 0
681	END DO
682
683
684	! Choose an integration bound well above wake top
685	! -----------------------------------------------------------------
686
687	! Determine Wake top pressure (Ptop) from buoyancy integral
688	! --------------------------------------------------------
689
690	Do i=1, klon
691	wk_adv(i) = .True.
692	Enddo
693	Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
694	dth, hw0, rho, delta_t_min, &
695	ktop, wk_adv, h_zzz, ptop1, ktop1)
696
697	!!print'("pkupper APPEL ",7i6)',0,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
698
699	IF (prt_level>=10) THEN
700	PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
701	ENDIF
702
703	! -5/ Set deltatw & deltaqw to 0 above kupper
704
705	DO k = 1, klev
706	DO i = 1, klon
707	IF (k>=kupper(i)) THEN
708	deltatw(i, k) = 0.
709	deltaqw(i, k) = 0.
710	d_deltatw2(i,k) = -deltatw0(i,k)
711	d_deltaqw2(i,k) = -deltaqw0(i,k)
712	END IF
713	END DO
714	END DO
715
716
717	! Vertical gradient of LS omega
718
719	DO k = 1, klev
720	DO i = 1, klon
721	IF (k<=kupper(i)) THEN
722	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
723	END IF
724	END DO
725	END DO
726
727	! Integrals (and wake top level number)
728	! --------------------------------------
729
730	! Initialize sum_thvx to 1st level virt. pot. temp.
731
732	DO i = 1, klon
733	z(i) = 1.
734	dz(i) = 1.
735	sum_thvx(i) = thx(i, 1)(1.+epsim1qx(i,1))*dz(i)
736	sum_dth(i) = 0.
737	END DO
738
739	DO k = 1, klev
740	DO i = 1, klon
741	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
742	IF (dz(i)>0) THEN
743	! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
744	z(i) = z(i) + dz(i)
745	sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
746	sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
747	sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
748	sum_thvx(i) = sum_thvx(i) + thx(i, k)(1.+epsim1qx(i,k))*dz(i)
749	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
750	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
751	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
752	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
753	END IF
754	END DO
755	END DO
756
757	DO i = 1, klon
758	hw0(i) = z(i)
759	END DO
760
761
762	! 2.1 - WAPE and mean forcing computation
763	! ---------------------------------------
764
765	! ---------------------------------------
766
767	! Means
768
769	DO i = 1, klon
770	av_thx(i) = sum_thx(i)/hw0(i)
771	av_tx(i) = sum_tx(i)/hw0(i)
772	av_qx(i) = sum_qx(i)/hw0(i)
773	av_thvx(i) = sum_thvx(i)/hw0(i)
774	! av_thve = sum_thve/hw0
775	av_dth(i) = sum_dth(i)/hw0(i)
776	av_dq(i) = sum_dq(i)/hw0(i)
777	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
778	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
779
780	wape(i) = -RGhw0(i)(av_dth(i)+ &
781	epsim1(av_thx(i)av_dq(i)+av_dth(i)av_qx(i)+av_dth(i)av_dq(i)))/av_thvx(i)
782
783	END DO
784	IF (CPPKEY_IOPHYS_WK) THEN
785	IF (.not.phys_sub) CALL iophys_ecrit('wape_a',1,'wape_a','J/kg',wape)
786	END IF
787
788	! 2.2 Prognostic variable update
789	! ------------------------------
790
791	! Filter out bad wakes
792
793	DO k = 1, klev
794	DO i = 1, klon
795	IF (wape(i)<0.) THEN
796	deltatw(i, k) = 0.
797	deltaqw(i, k) = 0.
798	dth(i, k) = 0.
799	d_deltatw2(i,k) = -deltatw0(i,k)
800	d_deltaqw2(i,k) = -deltaqw0(i,k)
801	END IF
802	END DO
803	END DO
804
805	DO i = 1, klon
806	IF (wape(i)<0.) THEN
807	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
808	sigmaw_targ = max(sigmad, sigd_con(i))
809	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
810	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
811	sigmaw(i) = sigmaw_targ
812	ENDIF !! (wape(i)<0.)
813	ENDDO
814	!
815	IF (iflag_wk_pop_dyn == 3) THEN
816	DO i = 1, klon
817	IF (wape(i)<0.) THEN
818	sigmaw_targ = max(sigmad, sigd_con(i))
819	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
820	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
821	asigmaw(i) = sigmaw_targ
822	ENDIF !! (wape(i)<0.)
823	ENDDO
824	ENDIF !! (iflag_wk_pop_dyn == 3)
825
826	DO i = 1, klon
827	IF (wape(i)<0.) THEN
828	wape(i) = 0.
829	cstar(i) = 0.
830	hw(i) = hwmin
831	fip(i) = 0.
832	gwake(i) = .FALSE.
833	ELSE
834	hw(i) = hw0(i)
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((qb(i,1)-sigmaw(i)*deltaqw(i,1)), &
845	(qb(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((qb(i,k)-sigmaw(i)*deltaqw(i,k)), &
850	(qb(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	alpha_tot(i) = 1.
861	END DO
862
863	IF (prt_level>=10) THEN
864	PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
865	sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
866	ENDIF
867
868
869	! C -----------------------------------------------------------------
870	! Sub-time-stepping
871	! -----------------
872
873	! wk_nsub and dtimesub definitions moved to begining of routine.
874	!! wk_nsub = 10
875	!! dtimesub = dtime/wk_nsub
876
877
878	! ------------------------------------------------------------------------
879	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
880	! ------------------------------------------------------------------------
881	!
882	DO isubstep = 1, wk_nsub
883	!
884	! ------------------------------------------------------------------------
885	! wk_adv is the logical flag enabling wake evolution in the time advance
886	! loop
887	DO i = 1, klon
888	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
889	END DO
890	IF (prt_level>=10) THEN
891	PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
892	isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
893
894	ENDIF
895
896	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
897	! negatif de ktop a kupper --------
898	! cc On calcule pour cela une densite wdens0 pour laquelle on
899	! aurait un entrainement nul ---
900	!jyg<
901	! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
902	! les poches sont insuffisantes pour accueillir tout le flux de masse
903	! des descentes unsaturees. Nous faisons alors l'hypothese que la
904	! convection profonde cree directement de nouvelles poches, sans passer
905	! par les thermiques. La nouvelle valeur de wdens est alors imposee.
906
907	DO i = 1, klon
908	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
909	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
910	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
911	IF ( iflag_wk_profile == 0 ) THEN
912	omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
913	RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
914	ELSE
915	omg(i, kupper(i)+1)=0.
916	ENDIF
917	wdens0 = (sigmaw(i)/(4.3.14)) &
918	((1.-sigmaw(i))omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))cstar(i)))**(2)
919	IF (prt_level >= 10) THEN
920	print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
921	omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
922	ENDIF
923	IF (wdens(i)<=wdens0*1.1) THEN
924	IF (iflag_wk_pop_dyn >= 1) THEN
925	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
926	d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
927	ENDIF
928	wdens(i) = wdens0
929	END IF
930	END IF
931	END DO
932
933	IF (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl) THEN
934	!!--------------------------------------------------------
935	!!Bug : computing gfl and rad_wk before changing sigmaw
936	!! This bug exists only for iflag_wk_pop_dyn=0. Otherwise, gfl and rad_wk
937	!! are computed within wake_popdyn
938	!!--------------------------------------------------------
939	DO i = 1, klon
940	IF (wk_adv(i)) THEN
941	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
942	rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
943	END IF
944	END DO
945	ENDIF ! (iflag_wk_pop_dyn == 0 .AND. ok_bug_gfl)
946	!!--------------------------------------------------------
947
948	DO i = 1, klon
949	IF (wk_adv(i)) THEN
950	sigmaw_targ = min(sigmaw(i), sigmaw_max)
951	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
952	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
953	sigmaw(i) = sigmaw_targ
954	END IF
955	END DO
956
957	IF (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl) THEN
958	!!--------------------------------------------------------
959	!!Fix : computing gfl and rad_wk after changing sigmaw
960	!!--------------------------------------------------------
961	DO i = 1, klon
962	IF (wk_adv(i)) THEN
963	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
964	rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
965	END IF
966	END DO
967	ENDIF ! (iflag_wk_pop_dyn == 0 .AND. .NOT.ok_bug_gfl)
968	!!--------------------------------------------------------
969
970	IF (iflag_wk_pop_dyn >= 1) THEN
971	! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
972	! Here, it has to be set to zero.
973	death_rate(:) = 0.
974	ENDIF
975
976	IF (iflag_wk_pop_dyn >= 3) THEN
977	DO i = 1, klon
978	IF (wk_adv(i)) THEN
979	sigmaw_targ = min(asigmaw(i), sigmaw_max)
980	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
981	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
982	asigmaw(i) = sigmaw_targ
983	ENDIF
984	ENDDO
985	ENDIF
986
987	!!--------------------------------------------------------
988	!!--------------------------------------------------------
989	IF (iflag_wk_pop_dyn == 1) THEN
990	!
991	CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
992	wdensmin, &
993	dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
994	d_awdens, d_wdens, d_sigmaw, &
995	iflag_wk_act, wk_adv, cin, wape, &
996	drdt, &
997	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
998	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
999	d_wdens_targ, d_sigmaw_targ)
1000
1001
1002	!!--------------------------------------------------------
1003	ELSEIF (iflag_wk_pop_dyn == 2) THEN
1004	!
1005	CALL wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
1006	wdensmin, &
1007	sigmaw, wdens, awdens, & !! state variables
1008	gfl, cstar, cin, wape, rad_wk, &
1009	d_sigmaw, d_wdens, d_awdens, & !! tendencies
1010	cont_fact, &
1011	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
1012	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
1013	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
1014	sigmaw=sigmaw-d_sigmaw
1015	wdens=wdens-d_wdens
1016	awdens=awdens-d_awdens
1017
1018	!!--------------------------------------------------------
1019	ELSEIF (iflag_wk_pop_dyn == 3) THEN
1020	IF (CPPKEY_IOPHYS_WK) THEN
1021	IF (phys_sub) THEN
1022	CALL iophys_ecrit('ptop',1,'ptop','Pa',ptop)
1023	CALL iophys_ecrit('sigmaw',1,'sigmaw','',sigmaw)
1024	CALL iophys_ecrit('asigmaw',1,'asigmaw','',asigmaw)
1025	CALL iophys_ecrit('wdens',1,'wdens','1/m2',wdens)
1026	CALL iophys_ecrit('awdens',1,'awdens','1/m2',awdens)
1027	CALL iophys_ecrit('rad_wk',1,'rad_wk','m',rad_wk)
1028	CALL iophys_ecrit('arad_wk',1,'arad_wk','m',arad_wk)
1029	CALL iophys_ecrit('irad_wk',1,'irad_wk','m',irad_wk)
1030	ENDIF
1031	END IF
1032	!
1033	CALL wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
1034	wdensmin, &
1035	sigmaw, asigmaw, wdens, awdens, & !! state variables
1036	gfl, agfl, cstar, cin, wape, &
1037	rad_wk, arad_wk, irad_wk, &
1038	d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies
1039	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
1040	d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
1041	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
1042	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
1043	sigmaw=sigmaw-d_sigmaw
1044	asigmaw=asigmaw-d_asigmaw
1045	wdens=wdens-d_wdens
1046	awdens=awdens-d_awdens
1047
1048	!!--------------------------------------------------------
1049	ELSEIF (iflag_wk_pop_dyn == 0) THEN
1050
1051	! cc nrlmd
1052
1053	DO i = 1, klon
1054	IF (wk_adv(i)) THEN
1055
1056	! cc nrlmd Introduction du taux de mortalite des poches et
1057	! test sur sigmaw_max=0.4
1058	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
1059	IF (sigmaw(i)>=sigmaw_max) THEN
1060	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
1061	ELSE
1062	death_rate(i) = 0.
1063	END IF
1064
1065	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
1066	dtimesub
1067	! $ - nat_rate(i)sigmaw(i)dtimesub
1068	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1069	! c $ death_rate(i),ktop(i),kupper(i)',
1070	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1071	! c $ death_rate(i),ktop(i),kupper(i)
1072
1073	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
1074	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
1075	! wdens = wdens0/(10.*sigmaw)
1076	! sigmaw =max(sigmaw,sigd_con)
1077	! sigmaw =max(sigmaw,sigmad)
1078	END IF
1079	END DO
1080
1081	ENDIF ! (iflag_wk_pop_dyn == 1) ... ELSEIF (iflag_wk_pop_dyn == 0)
1082	!!--------------------------------------------------------
1083	!!--------------------------------------------------------
1084
1085	IF (CPPKEY_IOPHYS_WK) THEN
1086	IF (phys_sub) THEN
1087	CALL iophys_ecrit('wdensa',1,'wdensa','m',wdens)
1088	CALL iophys_ecrit('awdensa',1,'awdensa','m',awdens)
1089	CALL iophys_ecrit('sigmawa',1,'sigmawa','m',sigmaw)
1090	CALL iophys_ecrit('asigmawa',1,'asigmawa','m',asigmaw)
1091	ENDIF
1092	END IF
1093	! calcul de la difference de vitesse verticale poche - zone non perturbee
1094	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
1095	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
1096	! IM 060208 au niveau k=1...
1097	!JYG 161013 Correction : maintenant omg est dimensionne a klev.
1098	DO k = 1, klev
1099	DO i = 1, klon
1100	IF (wk_adv(i)) THEN !!! nrlmd
1101	dp_deltomg(i, k) = 0.
1102	END IF
1103	END DO
1104	END DO
1105	DO k = 1, klev
1106	DO i = 1, klon
1107	IF (wk_adv(i)) THEN !!! nrlmd
1108	omg(i, k) = 0.
1109	END IF
1110	END DO
1111	END DO
1112
1113	DO i = 1, klon
1114	IF (wk_adv(i)) THEN
1115	z(i) = 0.
1116	omg(i, 1) = 0.
1117	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
1118	END IF
1119	END DO
1120
1121	DO k = 2, klev
1122	DO i = 1, klon
1123	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1124	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
1125	z(i) = z(i) + dz(i)
1126	dp_deltomg(i, k) = dp_deltomg(i, 1)
1127	omg(i, k) = dp_deltomg(i, 1)*z(i)
1128	END IF
1129	END DO
1130	END DO
1131
1132	DO i = 1, klon
1133	IF (wk_adv(i)) THEN
1134	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
1135	ztop(i) = z(i) + dztop(i)
1136	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
1137	END IF
1138	END DO
1139
1140	IF (prt_level>=10) THEN
1141	PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
1142	PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
1143	omgtop(igout), ptop(igout), ktop(igout)
1144	ENDIF
1145
1146	! -----------------
1147	! From m/s to Pa/s
1148	! -----------------
1149
1150	DO i = 1, klon
1151	IF (wk_adv(i)) THEN
1152	omgtop(i) = -rho(i, ktop(i))RGomgtop(i)
1153	!! LJYF dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
1154	dp_deltomg(i, 1) = omgtop(i)/min(ptop(i)-ph(i,1),-smallestreal)
1155	END IF
1156	END DO
1157
1158	DO k = 1, klev
1159	DO i = 1, klon
1160	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1161	omg(i, k) = -rho(i, k)RGomg(i, k)
1162	dp_deltomg(i, k) = dp_deltomg(i, 1)
1163	END IF
1164	END DO
1165	END DO
1166
1167	! raccordement lineaire de omg de ptop a pupper
1168
1169	DO i = 1, klon
1170	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
1171	IF ( iflag_wk_profile == 0 ) THEN
1172	omg(i, kupper(i)+1) =-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
1173	RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
1174	ELSE
1175	omg(i, kupper(i)+1) = 0.
1176	ENDIF
1177	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
1178	(ptop(i)-pupper(i))
1179	END IF
1180	END DO
1181
1182	! c DO i=1,klon
1183	! c print*,'Pente entre 0 et kupper (reference)'
1184	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
1185	! c print*,'Pente entre ktop et kupper'
1186	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
1187	! c ENDDO
1188	! c
1189	DO k = 1, klev
1190	DO i = 1, klon
1191	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
1192	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
1193	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
1194	END IF
1195	END DO
1196	END DO
1197	!! print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
1198	! cc nrlmd
1199	! c DO i=1,klon
1200	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
1201	! c END DO
1202	! cc
1203
1204
1205	! -- Compute wake average vertical velocity omgbw
1206
1207
1208	DO k = 1, klev
1209	DO i = 1, klon
1210	IF (wk_adv(i)) THEN
1211	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
1212	END IF
1213	END DO
1214	END DO
1215	! -- and its vertical gradient dp_omgbw
1216
1217	DO k = 1, klev-1
1218	DO i = 1, klon
1219	IF (wk_adv(i)) THEN
1220	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
1221	END IF
1222	END DO
1223	END DO
1224	DO i = 1, klon
1225	IF (wk_adv(i)) THEN
1226	dp_omgbw(i, klev) = 0.
1227	END IF
1228	END DO
1229
1230	! -- Upstream coefficients for omgb velocity
1231	! -- (alpha_up(k) is the coefficient of the value at level k)
1232	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
1233	DO k = 1, klev
1234	DO i = 1, klon
1235	IF (wk_adv(i)) THEN
1236	alpha_up(i, k) = 0.
1237	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
1238	END IF
1239	END DO
1240	END DO
1241
1242	! Matrix expressing [The,deltatw] from [Th1,Th2]
1243
1244	DO i = 1, klon
1245	IF (wk_adv(i)) THEN
1246	rre1(i) = 1. - sigmaw(i)
1247	rre2(i) = sigmaw(i)
1248	END IF
1249	END DO
1250	rrd1 = -1.
1251	rrd2 = 1.
1252
1253	! -- Get [Th1,Th2], dth and [q1,q2]
1254
1255	DO k = 1, klev
1256	DO i = 1, klon
1257	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1258	dth(i, k) = deltatw(i, k)/ppi(i, k)
1259	th1(i, k) = thb(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
1260	th2(i, k) = thb(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
1261	q1(i, k) = qb(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
1262	q2(i, k) = qb(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
1263	END IF
1264	END DO
1265	END DO
1266
1267	DO i = 1, klon
1268	IF (wk_adv(i)) THEN !!! nrlmd
1269	d_th1(i, 1) = 0.
1270	d_th2(i, 1) = 0.
1271	d_dth(i, 1) = 0.
1272	d_q1(i, 1) = 0.
1273	d_q2(i, 1) = 0.
1274	d_dq(i, 1) = 0.
1275	END IF
1276	END DO
1277
1278	DO k = 2, klev
1279	DO i = 1, klon
1280	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1281	d_th1(i, k) = th1(i, k-1) - th1(i, k)
1282	d_th2(i, k) = th2(i, k-1) - th2(i, k)
1283	d_dth(i, k) = dth(i, k-1) - dth(i, k)
1284	d_q1(i, k) = q1(i, k-1) - q1(i, k)
1285	d_q2(i, k) = q2(i, k-1) - q2(i, k)
1286	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
1287	END IF
1288	END DO
1289	END DO
1290
1291	DO i = 1, klon
1292	IF (wk_adv(i)) THEN
1293	omgbdth(i, 1) = 0.
1294	omgbdq(i, 1) = 0.
1295	END IF
1296	END DO
1297
1298	DO k = 2, klev
1299	DO i = 1, klon
1300	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
1301	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
1302	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
1303	END IF
1304	END DO
1305	END DO
1306
1307	!! IF (prt_level>=10) THEN
1308	IF (prt_level>=10 .and. wk_adv(igout)) THEN
1309	PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
1310	PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
1311	PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
1312	PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
1313	ENDIF
1314
1315	! -----------------------------------------------------------------
1316	DO k = 1, klev-1
1317	DO i = 1, klon
1318	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1319	! -----------------------------------------------------------------
1320
1321	! Compute redistribution (advective) term
1322
1323	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1324	(rrd1omg(i,k)sigmaw(i)*d_th1(i,k) - &
1325	rrd2omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1)- &
1326	(1.-alpha_up(i,k))*omgbdth(i,k)- &
1327	alpha_up(i,k+1)omgbdth(i,k+1))ppi(i, k)
1328	! print*,'d_d,k_ptop_provis(i)eltatw=', k, d_deltatw(i,k)
1329
1330	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1331	(rrd1omg(i,k)sigmaw(i)*d_q1(i,k)- &
1332	rrd2omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1)- &
1333	(1.-alpha_up(i,k))*omgbdq(i,k)- &
1334	alpha_up(i,k+1)*omgbdq(i,k+1))
1335	! print*,'d_deltaqw=', k, d_deltaqw(i,k)
1336
1337	! and increment large scale tendencies
1338
1339
1340
1341
1342	! C
1343	! -----------------------------------------------------------------
1344	d_tb(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)- &
1345	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/ &
1346	(ph(i,k)-ph(i,k+1)) &
1347	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
1348	(ph(i,k)-ph(i,k+1)) )*ppi(i, k)
1349
1350	d_qb(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)- &
1351	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/ &
1352	(ph(i,k)-ph(i,k+1)) &
1353	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
1354	(ph(i,k)-ph(i,k+1)) )
1355	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1356	d_tb(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
1357
1358	d_qb(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
1359
1360	END IF
1361	! cc
1362	END DO
1363	END DO
1364	! ------------------------------------------------------------------
1365
1366	IF (prt_level>=10) THEN
1367	PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
1368	PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
1369	ENDIF
1370
1371	! Increment state variables
1372	!jyg<
1373	IF (iflag_wk_pop_dyn >= 1) THEN
1374	DO k = 1, klev
1375	DO i = 1, klon
1376	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1377	detr(i,k) = - d_sig_death(i) - d_sig_col(i)
1378	entr(i,k) = d_sig_gen(i)
1379	ENDIF
1380	ENDDO
1381	ENDDO
1382	ELSE ! (iflag_wk_pop_dyn >= 1)
1383	DO k = 1, klev
1384	DO i = 1, klon
1385	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1386	detr(i, k) = 0.
1387
1388	entr(i, k) = 0.
1389	ENDIF
1390	ENDDO
1391	ENDDO
1392	ENDIF ! (iflag_wk_pop_dyn >= 1)
1393
1394
1395
1396	DO k = 1, klev
1397	DO i = 1, klon
1398	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1399	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1400	! cc
1401
1402
1403
1404	! Coefficient de repartition
1405
1406	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1407	(ph(i,kupper(i))-ph(i,1))
1408	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
1409	(ph(i,1)-ph(i,kupper(i)))
1410
1411
1412	! Reintroduce compensating subsidence term.
1413
1414	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1415	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1416	! . /(1-sigmaw)
1417	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1418	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1419	! . /(1-sigmaw)
1420
1421	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
1422	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
1423	! . /(1-sigmaw)
1424	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
1425	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
1426	! . /(1-sigmaw)
1427
1428	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1429	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1430	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1431
1432	!
1433
1434	! cc nrlmd Prise en compte du taux de mortalite
1435	! cc Definitions de entr, detr
1436	!jyg<
1437	!! detr(i, k) = 0.
1438	!!
1439	!! entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1440	!! sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1441	!!
1442	entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
1443	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
1444	!>jyg
1445	wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1446
1447	! cc wkspread(i,k) =
1448	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
1449	! cc $ sigmaw(i)
1450
1451
1452	! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
1453	! Jingmei
1454
1455	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1456	! & Tgw(i,k),deltatw(i,k)
1457	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
1458	dtimesub
1459	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1460	ff(i) = d_deltatw(i, k)/dtimesub
1461
1462	! Sans GW
1463
1464	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-wkspread(k)deltatw(k))
1465
1466	! GW formule 1
1467
1468	! deltatw(k) = deltatw(k)+dtimesub*
1469	! $ (ff+dtKE(k) - wkspread(k)deltatw(k)-Tgw(k)deltatw(k))
1470
1471	! GW formule 2
1472
1473	IF (dtimesub*tgw(i,k)<1.E-10) THEN
1474	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
1475	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1476	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
1477	tgw(i,k)*deltatw(i,k) )
1478	ELSE
1479	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i,k)))* &
1480	(ff(i)+dtke(i,k) - &
1481	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1482	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - &
1483	tgw(i,k)*deltatw(i,k) )
1484	END IF
1485
1486	dth(i, k) = deltatw(i, k)/ppi(i, k)
1487
1488	gg(i) = d_deltaqw(i, k)/dtimesub
1489
1490	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
1491	entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
1492	(death_rate(i)sigmaw(i)+detr(i,k))deltaqw(i,k)/(1.-sigmaw(i)))
1493	! cc
1494
1495	! cc nrlmd
1496	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1497	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1498	! cc
1499	END IF
1500	END DO
1501	END DO
1502
1503
1504	! Scale tendencies so that water vapour remains positive in w and x.
1505
1506	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
1507	d_deltaqw, sigmaw, d_sigmaw, alpha)
1508	!
1509	! Alpha_tot = Product of all the alpha's
1510	DO i = 1, klon
1511	IF (wk_adv(i)) THEN
1512	alpha_tot(i) = alpha_tot(i)*alpha(i)
1513	END IF
1514	END DO
1515
1516	! cc nrlmd
1517	! c print*,'alpha'
1518	! c do i=1,klon
1519	! c print*,alpha(i)
1520	! c end do
1521	! cc
1522	DO k = 1, klev
1523	DO i = 1, klon
1524	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1525	d_tb(i, k) = alpha(i)*d_tb(i, k)
1526	d_qb(i, k) = alpha(i)*d_qb(i, k)
1527	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1528	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1529	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1530	END IF
1531	END DO
1532	END DO
1533	DO i = 1, klon
1534	IF (wk_adv(i)) THEN
1535	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1536	END IF
1537	END DO
1538
1539	! Update large scale variables and wake variables
1540	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1541	! IM 060208 DO k = 1,kupper(i)
1542	DO k = 1, klev
1543	DO i = 1, klon
1544	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1545	dtls(i, k) = dtls(i, k) + d_tb(i, k)
1546	dqls(i, k) = dqls(i, k) + d_qb(i, k)
1547	! cc nrlmd
1548	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1549	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1550	! cc
1551	END IF
1552	END DO
1553	END DO
1554	DO k = 1, klev
1555	DO i = 1, klon
1556	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1557	tb(i, k) = tb0(i, k) + dtls(i, k)
1558	qb(i, k) = qb0(i, k) + dqls(i, k)
1559	thb(i, k) = tb(i, k)/ppi(i, k)
1560	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1561	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1562	dth(i, k) = deltatw(i, k)/ppi(i, k)
1563	! c print,'k,qx,qw',k,qb(i,k)-sigmaw(i)deltaqw(i,k)
1564	! c $ ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1565	END IF
1566	END DO
1567	END DO
1568	!
1569	DO i = 1, klon
1570	IF (wk_adv(i)) THEN
1571	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1572	d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
1573	END IF
1574	END DO
1575	!jyg<
1576	IF (iflag_wk_pop_dyn >= 1) THEN
1577	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! sigmaw !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1578	! Cumulatives
1579	DO i = 1, klon
1580	IF (wk_adv(i)) THEN
1581	d_sig_gen2(i) = d_sig_gen2(i) + d_sig_gen(i)
1582	d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
1583	d_sig_col2(i) = d_sig_col2(i) + d_sig_col(i)
1584	d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
1585	d_sig_bnd2(i) = d_sig_bnd2(i) + d_sig_bnd(i)
1586	END IF
1587	END DO
1588	! Bounds
1589	DO i = 1, klon
1590	IF (wk_adv(i)) THEN
1591	sigmaw_targ = max(sigmaw(i),sigmad)
1592	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
1593	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1594	sigmaw(i) = sigmaw_targ
1595	END IF
1596	END DO
1597	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! wdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1598	! Cumulatives
1599	DO i = 1, klon
1600	IF (wk_adv(i)) THEN
1601	wdens(i) = wdens(i) + d_wdens(i)
1602	d_wdens2(i) = d_wdens2(i) + d_wdens(i)
1603	d_dens_gen2(i) = d_dens_gen2(i) + d_dens_gen(i)
1604	d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
1605	d_dens_col2(i) = d_dens_col2(i) + d_dens_col(i)
1606	d_dens_bnd2(i) = d_dens_bnd2(i) + d_dens_bnd(i)
1607	END IF
1608	END DO
1609	! Bounds
1610	DO i = 1, klon
1611	IF (wk_adv(i)) THEN
1612	wdens_targ = max(wdens(i),wdensmin)
1613	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
1614	d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
1615	wdens(i) = wdens_targ
1616	END IF
1617	END DO
1618	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1619	! Cumulatives
1620	DO i = 1, klon
1621	IF (wk_adv(i)) THEN
1622	awdens(i) = awdens(i) + d_awdens(i)
1623	d_awdens2(i) = d_awdens2(i) + d_awdens(i)
1624	END IF
1625	END DO
1626	! Bounds
1627	DO i = 1, klon
1628	IF (wk_adv(i)) THEN
1629	wdens_targ = min( max(awdens(i),0.), wdens(i) )
1630	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
1631	d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
1632	awdens(i) = wdens_targ
1633	END IF
1634	END DO
1635	!
1636	IF (iflag_wk_pop_dyn >= 2) THEN
1637	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! awdens again for iflag_wk_pop_dyn >= 2!!!!!!
1638	! Cumulatives
1639	DO i = 1, klon
1640	IF (wk_adv(i)) THEN
1641	d_adens_death2(i) = d_adens_death2(i) + d_adens_death(i)
1642	d_adens_icol2(i) = d_adens_icol2(i) + d_adens_icol(i)
1643	d_adens_acol2(i) = d_adens_acol2(i) + d_adens_acol(i)
1644	d_adens_bnd2(i) = d_adens_bnd2(i) + d_adens_bnd(i)
1645	END IF
1646	END DO
1647	! Bounds
1648	DO i = 1, klon
1649	IF (wk_adv(i)) THEN
1650	wdens_targ = min( max(awdens(i),0.), wdens(i) )
1651	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
1652	awdens(i) = wdens_targ
1653	END IF
1654	END DO
1655	!
1656	IF (iflag_wk_pop_dyn == 3) THEN
1657	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! asigmaw for iflag_wk_pop_dyn = 3!!!!!!
1658	! Cumulatives
1659	DO i = 1, klon
1660	IF (wk_adv(i)) THEN
1661	asigmaw(i) = asigmaw(i) + d_asigmaw(i)
1662	d_asigmaw2(i) = d_asigmaw2(i) + d_asigmaw(i)
1663	d_asig_death2(i) = d_asig_death2(i) + d_asig_death(i)
1664	d_asig_spread2(i) = d_asig_spread2(i) + d_asig_spread(i)
1665	d_asig_iicol2(i) = d_asig_iicol2(i) + d_asig_iicol(i)
1666	d_asig_aicol2(i) = d_asig_aicol2(i) + d_asig_aicol(i)
1667	d_asig_bnd2(i) = d_asig_bnd2(i) + d_asig_bnd(i)
1668	END IF
1669	END DO
1670	! Bounds
1671	DO i = 1, klon
1672	IF (wk_adv(i)) THEN
1673	! asigmaw lower bound set to sigmad/2 in order to allow asigmaw values lower than sigmad.
1674	!! sigmaw_targ = min(max(asigmaw(i),sigmad),sigmaw(i))
1675	sigmaw_targ = min(max(asigmaw(i),sigmad/2.),sigmaw(i))
1676	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
1677	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
1678	asigmaw(i) = sigmaw_targ
1679	END IF
1680	END DO
1681
1682	IF (CPPKEY_IOPHYS_WK) THEN
1683	IF (phys_sub) THEN
1684	CALL iophys_ecrit('wdensb',1,'wdensb','m',wdens)
1685	CALL iophys_ecrit('awdensb',1,'awdensb','m',awdens)
1686	CALL iophys_ecrit('sigmawb',1,'sigmawb','m',sigmaw)
1687	CALL iophys_ecrit('asigmawb',1,'asigmawb','m',asigmaw)
1688	!
1689	call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
1690	call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
1691	call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
1692	call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
1693	call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
1694	!
1695	call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
1696	call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
1697	call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
1698	call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
1699	call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
1700	!
1701	CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
1702	CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
1703	CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
1704	CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
1705	CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
1706	CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
1707	!
1708	CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
1709	CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
1710	CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
1711	CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
1712	CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
1713	CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
1714	ENDIF
1715	END IF
1716	ENDIF ! (iflag_wk_pop_dyn == 3)
1717	ENDIF ! (iflag_wk_pop_dyn >= 2)
1718	ENDIF ! (iflag_wk_pop_dyn >= 1)
1719
1720
1721
1722	Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
1723	dth, hw, rho, delta_t_min, &
1724	ktop, wk_adv, h_zzz, ptop1, ktop1)
1725	!! print'("pkupper APPEL ",7i6)',isubstep,int(ptop/100.),int(ptop1/100.),int(pupper/100.),ktop,ktop1,kupper
1726
1727	! 5/ Set deltatw & deltaqw to 0 above kupper
1728
1729	DO k = 1, klev
1730	DO i = 1, klon
1731	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
1732	deltatw(i, k) = 0.
1733	deltaqw(i, k) = 0.
1734	d_deltatw2(i,k) = -deltatw0(i,k)
1735	d_deltaqw2(i,k) = -deltaqw0(i,k)
1736	END IF
1737	END DO
1738	END DO
1739
1740
1741	! -------------Cstar computation---------------------------------
1742	DO i = 1, klon
1743	IF (wk_adv(i)) THEN !!! nrlmd
1744	sum_thx(i) = 0.
1745	sum_tx(i) = 0.
1746	sum_qx(i) = 0.
1747	sum_thvx(i) = 0.
1748	sum_dth(i) = 0.
1749	sum_dq(i) = 0.
1750	sum_dtdwn(i) = 0.
1751	sum_dqdwn(i) = 0.
1752
1753	av_thx(i) = 0.
1754	av_tx(i) = 0.
1755	av_qx(i) = 0.
1756	av_thvx(i) = 0.
1757	av_dth(i) = 0.
1758	av_dq(i) = 0.
1759	av_dtdwn(i) = 0.
1760	av_dqdwn(i) = 0.
1761	END IF
1762	END DO
1763
1764	! Integrals (and wake top level number)
1765	! --------------------------------------
1766
1767	! Initialize sum_thvx to 1st level virt. pot. temp.
1768
1769	DO i = 1, klon
1770	IF (wk_adv(i)) THEN !!! nrlmd
1771	z(i) = 1.
1772	dz(i) = 1.
1773	sum_thvx(i) = thx(i, 1)(1.+epsim1qx(i,1))*dz(i)
1774	sum_dth(i) = 0.
1775	END IF
1776	END DO
1777
1778	DO k = 1, klev
1779	DO i = 1, klon
1780	IF (wk_adv(i)) THEN !!! nrlmd
1781	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
1782	IF (dz(i)>0) THEN
1783	z(i) = z(i) + dz(i)
1784	sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
1785	sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
1786	sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
1787	sum_thvx(i) = sum_thvx(i) + thx(i, k)(1.+epsim1qx(i,k))*dz(i)
1788	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1789	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1790	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1791	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1792	END IF
1793	END IF
1794	END DO
1795	END DO
1796
1797	DO i = 1, klon
1798	IF (wk_adv(i)) THEN !!! nrlmd
1799	hw0(i) = z(i)
1800	END IF
1801	END DO
1802
1803
1804	! - WAPE and mean forcing computation
1805	! ---------------------------------------
1806
1807	! ---------------------------------------
1808
1809	! Means
1810
1811	DO i = 1, klon
1812	IF (wk_adv(i)) THEN !!! nrlmd
1813	av_thx(i) = sum_thx(i)/hw0(i)
1814	av_tx(i) = sum_tx(i)/hw0(i)
1815	av_qx(i) = sum_qx(i)/hw0(i)
1816	av_thvx(i) = sum_thvx(i)/hw0(i)
1817	av_dth(i) = sum_dth(i)/hw0(i)
1818	av_dq(i) = sum_dq(i)/hw0(i)
1819	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1820	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1821
1822	wape(i) = -RGhw0(i)(av_dth(i)+epsim1(av_thx(i)av_dq(i) + &
1823	av_dth(i)av_qx(i)+av_dth(i)av_dq(i)))/av_thvx(i)
1824	END IF
1825	END DO
1826
1827
1828	! Filter out bad wakes
1829
1830	DO k = 1, klev
1831	DO i = 1, klon
1832	IF (wk_adv(i)) THEN !!! nrlmd
1833	IF (wape(i)<0.) THEN
1834	deltatw(i, k) = 0.
1835	deltaqw(i, k) = 0.
1836	dth(i, k) = 0.
1837	d_deltatw2(i,k) = -deltatw0(i,k)
1838	d_deltaqw2(i,k) = -deltaqw0(i,k)
1839	END IF
1840	END IF
1841	END DO
1842	END DO
1843
1844	DO i = 1, klon
1845	IF (wk_adv(i)) THEN !!! nrlmd
1846	IF (wape(i)<0.) THEN
1847	wape(i) = 0.
1848	cstar(i) = 0.
1849	hw(i) = hwmin
1850	!jyg<
1851	!! sigmaw(i) = max(sigmad, sigd_con(i))
1852	sigmaw_targ = max(sigmad, sigd_con(i))
1853	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
1854	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1855	sigmaw(i) = sigmaw_targ
1856	!
1857	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
1858	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
1859	asigmaw(i) = sigmaw_targ
1860	!>jyg
1861	fip(i) = 0.
1862	gwake(i) = .FALSE.
1863	ELSE
1864	cstar(i) = starksqrt(2.wape(i))
1865	gwake(i) = .TRUE.
1866	END IF
1867	END IF
1868	END DO
1869	!
1870	! ------------------------------------------------------------------------
1871	!
1872	END DO ! isubstep end sub-timestep loop
1873	!
1874	! ------------------------------------------------------------------------
1875	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1876	! ------------------------------------------------------------------------
1877	!
1878
1879	IF (CPPKEY_IOPHYS_WK) THEN
1880	IF (.not.phys_sub) CALL iophys_ecrit('wape_b',1,'wape_b','J/kg',wape)
1881	END IF
1882	IF (prt_level>=10) THEN
1883	PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
1884	sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
1885	ENDIF
1886
1887
1888	! ----------------------------------------------------------
1889	! Determine wake final state; recompute wape, cstar, ktop;
1890	! filter out bad wakes.
1891	! ----------------------------------------------------------
1892
1893	! 2.1 - Undisturbed area and Wake integrals
1894	! ---------------------------------------------------------
1895
1896	DO i = 1, klon
1897	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
1898	IF (ok_qx_qw(i)) THEN
1899	! cc
1900	z(i) = 0.
1901	sum_thx(i) = 0.
1902	sum_tx(i) = 0.
1903	sum_qx(i) = 0.
1904	sum_thvx(i) = 0.
1905	sum_dth(i) = 0.
1906	sum_half_dth(i) = 0.
1907	sum_dq(i) = 0.
1908	sum_dtdwn(i) = 0.
1909	sum_dqdwn(i) = 0.
1910
1911	av_thx(i) = 0.
1912	av_tx(i) = 0.
1913	av_qx(i) = 0.
1914	av_thvx(i) = 0.
1915	av_dth(i) = 0.
1916	av_dq(i) = 0.
1917	av_dtdwn(i) = 0.
1918	av_dqdwn(i) = 0.
1919
1920	dthmin(i) = -delta_t_min
1921	END IF
1922	END DO
1923	! Potential temperatures and humidity
1924	! ----------------------------------------------------------
1925
1926	DO k = 1, klev
1927	DO i = 1, klon
1928	! cc nrlmd IF ( wk_adv(i)) THEN
1929	IF (ok_qx_qw(i)) THEN
1930	! cc
1931	rho(i, k) = p(i, k)/(RD*tb(i,k))
1932	IF (k==1) THEN
1933	rhoh(i, k) = ph(i, k)/(RD*tb(i,k))
1934	zhh(i, k) = 0
1935	ELSE
1936	rhoh(i, k) = ph(i, k)2./(RD(tb(i,k)+tb(i,k-1)))
1937	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
1938	END IF
1939	thb(i, k) = tb(i, k)/ppi(i, k)
1940	thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
1941	tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i)
1942	qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i)
1943	dth(i, k) = deltatw(i, k)/ppi(i, k)
1944	END IF
1945	END DO
1946	END DO
1947
1948	! Integrals (and wake top level number)
1949	! -----------------------------------------------------------
1950
1951	! Initialize sum_thvx to 1st level virt. pot. temp.
1952
1953	DO i = 1, klon
1954	! cc nrlmd IF ( wk_adv(i)) THEN
1955	IF (ok_qx_qw(i)) THEN
1956	! cc
1957	z(i) = 1.
1958	dz(i) = 1.
1959	dz_half(i) = 1.
1960	sum_thvx(i) = thx(i, 1)(1.+epsim1qx(i,1))*dz(i)
1961	sum_dth(i) = 0.
1962	END IF
1963	END DO
1964
1965	DO k = 1, klev
1966	DO i = 1, klon
1967	! cc nrlmd IF ( wk_adv(i)) THEN
1968	IF (ok_qx_qw(i)) THEN
1969	! cc
1970	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
1971	dz_half(i) = -(amax1(ph(i,k+1),0.5(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)RG)
1972	IF (dz(i)>0) THEN
1973	z(i) = z(i) + dz(i)
1974	sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i)
1975	sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i)
1976	sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i)
1977	sum_thvx(i) = sum_thvx(i) + thx(i, k)(1.+epsim1qx(i,k))*dz(i)
1978	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1979	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
1980	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
1981	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
1982	!
1983	dthmin(i) = min(dthmin(i), dth(i,k))
1984	END IF
1985	IF (dz_half(i)>0) THEN
1986	sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
1987	END IF
1988	END IF
1989	END DO
1990	END DO
1991
1992	DO i = 1, klon
1993	! cc nrlmd IF ( wk_adv(i)) THEN
1994	IF (ok_qx_qw(i)) THEN
1995	! cc
1996	hw0(i) = z(i)
1997	END IF
1998	END DO
1999
2000	! - WAPE and mean forcing computation
2001	! -------------------------------------------------------------
2002
2003	! Means
2004
2005	DO i = 1, klon
2006	! cc nrlmd IF ( wk_adv(i)) THEN
2007	IF (ok_qx_qw(i)) THEN
2008	! cc
2009	av_thx(i) = sum_thx(i)/hw0(i)
2010	av_tx(i) = sum_tx(i)/hw0(i)
2011	av_qx(i) = sum_qx(i)/hw0(i)
2012	av_thvx(i) = sum_thvx(i)/hw0(i)
2013	av_dth(i) = sum_dth(i)/hw0(i)
2014	av_dq(i) = sum_dq(i)/hw0(i)
2015	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
2016	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
2017
2018	wape2(i) = -RGhw0(i)(av_dth(i)+epsim1(av_thx(i)av_dq(i) + &
2019	av_dth(i)av_qx(i)+av_dth(i)av_dq(i)))/av_thvx(i)
2020	END IF
2021	END DO
2022	IF (CPPKEY_IOPHYS_WK) THEN
2023	IF (.not.phys_sub) CALL iophys_ecrit('wape2_a',1,'wape2_a','J/kg',wape2)
2024	END IF
2025
2026
2027	! Prognostic variable update
2028	! ------------------------------------------------------------
2029
2030	! Filter out bad wakes
2031
2032	IF (iflag_wk_check_trgl>=1) THEN
2033	! Check triangular shape of dth profile
2034	DO i = 1, klon
2035	IF (ok_qx_qw(i)) THEN
2036	!! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
2037	!! print ,'wake, 2.sum_dth(i)/(hw0(i)*dthmin(i)) ', &
2038	!! 2.sum_dth(i)/(hw0(i)dthmin(i))
2039	!! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
2040	!! sum_half_dth(i), sum_dth(i)
2041	IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
2042	wape2(i) = -1.
2043	!! print *,'wake, rej 1'
2044	ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.sum_dth(i)/(hw0(i)dthmin(i)) - 1.) > 0.5) THEN
2045	wape2(i) = -1.
2046	!! print *,'wake, rej 2'
2047	ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
2048	wape2(i) = -1.
2049	!! print *,'wake, rej 3'
2050	END IF
2051	END IF
2052	END DO
2053	END IF
2054	IF (CPPKEY_IOPHYS_WK) THEN
2055	IF (.not.phys_sub) CALL iophys_ecrit('wape2_b',1,'wape2_b','J/kg',wape2)
2056	END IF
2057
2058
2059	DO k = 1, klev
2060	DO i = 1, klon
2061	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
2062	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
2063	! cc
2064	deltatw(i, k) = 0.
2065	deltaqw(i, k) = 0.
2066	dth(i, k) = 0.
2067	d_deltatw2(i,k) = -deltatw0(i,k)
2068	d_deltaqw2(i,k) = -deltaqw0(i,k)
2069	END IF
2070	END DO
2071	END DO
2072
2073
2074	DO i = 1, klon
2075	! cc nrlmd IF ( wk_adv(i)) THEN
2076	IF (ok_qx_qw(i)) THEN
2077	! cc
2078	IF (wape2(i)<0.) THEN
2079	wape2(i) = 0.
2080	cstar2(i) = 0.
2081	hw(i) = hwmin
2082	!jyg<
2083	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
2084	sigmaw_targ = max(sigmad, sigd_con(i))
2085	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
2086	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2087	sigmaw(i) = sigmaw_targ
2088	!
2089	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
2090	d_asigmaw2(i) = d_asigmaw2(i) + sigmaw_targ - asigmaw(i)
2091	asigmaw(i) = sigmaw_targ
2092	!>jyg
2093	fip(i) = 0.
2094	gwake(i) = .FALSE.
2095	ELSE
2096	IF (prt_level>=10) PRINT *, 'wape2>0'
2097	cstar2(i) = starksqrt(2.wape2(i))
2098	gwake(i) = .TRUE.
2099	END IF
2100	IF (CPPKEY_IOPHYS_WK) THEN
2101	IF (.not.phys_sub) CALL iophys_ecrit('cstar2',1,'cstar2','J/kg',cstar2)
2102	END IF
2103	END IF ! (ok_qx_qw(i))
2104	END DO
2105
2106	DO i = 1, klon
2107	! cc nrlmd IF ( wk_adv(i)) THEN
2108	IF (ok_qx_qw(i)) THEN
2109	! cc
2110	ktopw(i) = ktop(i)
2111	END IF
2112	END DO
2113
2114	DO i = 1, klon
2115	! cc nrlmd IF ( wk_adv(i)) THEN
2116	IF (ok_qx_qw(i)) THEN
2117	! cc
2118	IF (ktopw(i)>0 .AND. gwake(i)) THEN
2119
2120	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
2121	! cc heff = 600.
2122	! Utilisation de la hauteur hw
2123	! c heff = 0.7*hw
2124	heff(i) = hw(i)
2125
2126	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
2127	sqrt(sigmaw(i)wdens(i)3.14)
2128	fip(i) = alpk*fip(i)
2129	! jyg2
2130	ELSE
2131	fip(i) = 0.
2132	END IF
2133	END IF
2134	END DO
2135	IF (iflag_wk_pop_dyn >= 3) THEN
2136	IF (CPPKEY_IOPHYS_WK) THEN
2137	IF (.not.phys_sub) THEN
2138	CALL iophys_ecrit('fip',1,'fip','J/kg',fip)
2139	CALL iophys_ecrit('hw',1,'hw','J/kg',hw)
2140	CALL iophys_ecrit('ptop',1,'ptop','J/kg',ptop)
2141	CALL iophys_ecrit('wdens',1,'wdens','J/kg',wdens)
2142	CALL iophys_ecrit('awdens',1,'awdens','m',awdens)
2143	CALL iophys_ecrit('sigmaw',1,'sigmaw','m',sigmaw)
2144	CALL iophys_ecrit('asigmaw',1,'asigmaw','m',asigmaw)
2145	!
2146	CALL iophys_ecrit('rad_wk',1,'rad_wk','J/kg',rad_wk)
2147	CALL iophys_ecrit('arad_wk',1,'arad_wk','J/kg',arad_wk)
2148	CALL iophys_ecrit('irad_wk',1,'irad_wk','J/kg',irad_wk)
2149	!
2150	call iophys_ecrit('d_wdens2',1,'d_wdens2','',d_wdens2)
2151	call iophys_ecrit('d_dens_gen2',1,'d_dens_gen2','',d_dens_gen2)
2152	call iophys_ecrit('d_dens_death2',1,'d_dens_death2','',d_dens_death2)
2153	call iophys_ecrit('d_dens_col2',1,'d_dens_col2','',d_dens_col2)
2154	call iophys_ecrit('d_dens_bnd2',1,'d_dens_bnd2','',d_dens_bnd2)
2155	!
2156	call iophys_ecrit('d_awdens2',1,'d_awdens2','',d_awdens2)
2157	call iophys_ecrit('d_adens_death2',1,'d_adens_death2','',d_adens_death2)
2158	call iophys_ecrit('d_adens_icol2',1,'d_adens_icol2','',d_adens_icol2)
2159	call iophys_ecrit('d_adens_acol2',1,'d_adens_acol2','',d_adens_acol2)
2160	call iophys_ecrit('d_adens_bnd2',1,'d_adens_bnd2','',d_adens_bnd2)
2161	!
2162	CALL iophys_ecrit('d_sigmaw2',1,'d_sigmaw2','',d_sigmaw2)
2163	CALL iophys_ecrit('d_sig_gen2',1,'d_sig_gen2','m',d_sig_gen2)
2164	CALL iophys_ecrit('d_sig_spread2',1,'d_sig_spread2','',d_sig_spread2)
2165	CALL iophys_ecrit('d_sig_col2',1,'d_sig_col2','',d_sig_col2)
2166	CALL iophys_ecrit('d_sig_death2',1,'d_sig_death2','',d_sig_death2)
2167	CALL iophys_ecrit('d_sig_bnd2',1,'d_sig_bnd2','',d_sig_bnd2)
2168	!
2169	CALL iophys_ecrit('d_asigmaw2',1,'d_asigmaw2','',d_asigmaw2)
2170	CALL iophys_ecrit('d_asig_spread2',1,'d_asig_spread2','m',d_asig_spread2)
2171	CALL iophys_ecrit('d_asig_aicol2',1,'d_asig_aicol2','m',d_asig_aicol2)
2172	CALL iophys_ecrit('d_asig_iicol2',1,'d_asig_iicol2','m',d_asig_iicol2)
2173	CALL iophys_ecrit('d_asig_death2',1,'d_asig_death2','m',d_asig_death2)
2174	CALL iophys_ecrit('d_asig_bnd2',1,'d_asig_bnd2','m',d_asig_bnd2)
2175	ENDIF ! (.not.phys_sub)
2176	END IF
2177	ENDIF ! (iflag_wk_pop_dyn >= 3)
2178	! Limitation de sigmaw
2179
2180	! cc nrlmd
2181	! DO i=1,klon
2182	! IF (OK_qx_qw(i)) THEN
2183	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
2184	! ENDIF
2185	! ENDDO
2186	! cc
2187
2188	!jyg<
2189	IF (iflag_wk_pop_dyn >= 1) THEN
2190	DO i = 1, klon
2191	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2192	.NOT. ok_qx_qw(i) .OR. (wdens(i) < wdensthreshold)
2193	!! .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
2194	ENDDO
2195	ELSE ! (iflag_wk_pop_dyn >= 1)
2196	DO i = 1, klon
2197	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2198	.NOT. ok_qx_qw(i)
2199	ENDDO
2200	ENDIF ! (iflag_wk_pop_dyn >= 1)
2201	!>jyg
2202
2203	DO k = 1, klev
2204	DO i = 1, klon
2205	!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2206	!!jyg .NOT. ok_qx_qw(i)) THEN
2207	IF (kill_wake(i)) THEN
2208	! cc
2209	dtls(i, k) = 0.
2210	dqls(i, k) = 0.
2211	deltatw(i, k) = 0.
2212	deltaqw(i, k) = 0.
2213	d_deltatw2(i,k) = -deltatw0(i,k)
2214	d_deltaqw2(i,k) = -deltaqw0(i,k)
2215	END IF ! (kill_wake(i))
2216	END DO
2217	END DO
2218
2219	DO i = 1, klon
2220	!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2221	!!jyg .NOT. ok_qx_qw(i)) THEN
2222	IF (kill_wake(i)) THEN
2223	ktopw(i) = 0
2224	wape(i) = 0.
2225	cstar(i) = 0.
2226	!!jyg Outside subroutine "Wake" hw, wdens sigmaw and asigmaw are zero when there are no wakes
2227	!! hw(i) = hwmin !jyg
2228	!! sigmaw(i) = sigmad !jyg
2229	hw(i) = 0. !jyg
2230	fip(i) = 0.
2231	!
2232	!! sigmaw(i) = 0. !jyg
2233	sigmaw_targ = 0.
2234	d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
2235	!! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2236	d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i) ! _in = correction jyg 20220124
2237	sigmaw(i) = sigmaw_targ
2238	!
2239	IF (iflag_wk_pop_dyn >= 3) THEN
2240	sigmaw_targ = 0.
2241	d_asig_bnd2(i) = d_asig_bnd2(i) + sigmaw_targ - asigmaw(i)
2242	!! d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2243	d_asigmaw2(i) = sigmaw_targ - asigmaw_in(i) ! _in = correction jyg 20220124
2244	asigmaw(i) = sigmaw_targ
2245	ELSE
2246	asigmaw(i) = 0.
2247	ENDIF ! (iflag_wk_pop_dyn >= 3)
2248	!
2249	IF (iflag_wk_pop_dyn >= 1) THEN
2250	!! awdens(i) = 0.
2251	!! wdens(i) = 0.
2252	wdens_targ = 0.
2253	d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
2254	!! d_wdens2(i) = wdens_targ - wdens(i)
2255	d_wdens2(i) = wdens_targ - wdens_in(i) ! jyg 20220916
2256	wdens(i) = wdens_targ
2257	wdens_targ = 0.
2258	!!jyg: bug fix : the d_adens_bnd2 computation must be before the update of awdens.
2259	IF (iflag_wk_pop_dyn >= 2) THEN
2260	d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
2261	ENDIF ! (iflag_wk_pop_dyn >= 2)
2262	!! d_awdens2(i) = wdens_targ - awdens(i)
2263	d_awdens2(i) = wdens_targ - awdens_in(i) ! jyg 20220916
2264	awdens(i) = wdens_targ
2265	!! IF (iflag_wk_pop_dyn == 2) THEN
2266	!! d_adens_bnd2(i) = d_adens_bnd2(i) + wdens_targ - awdens(i)
2267	!! ENDIF ! (iflag_wk_pop_dyn == 2)
2268	ENDIF ! (iflag_wk_pop_dyn >= 1)
2269	ELSE ! (kill_wake(i))
2270	wape(i) = wape2(i)
2271	cstar(i) = cstar2(i)
2272	END IF ! (kill_wake(i))
2273	! c print*,'wape wape2 ktopw OK_qx_qw =',
2274	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
2275	END DO
2276
2277	IF (prt_level>=10) THEN
2278	PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
2279	wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
2280	ENDIF
2281	IF (CPPKEY_IOPHYS_WK) THEN
2282	IF (.not.phys_sub) CALL iophys_ecrit('wape_c',1,'wape_c','J/kg',wape)
2283	END IF
2284
2285
2286	! -----------------------------------------------------------------
2287	! Get back to tendencies per second
2288
2289	DO k = 1, klev
2290	DO i = 1, klon
2291
2292	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
2293	!jyg<
2294	!! IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
2295	IF (ok_qx_qw(i)) THEN
2296	!>jyg
2297	! cc
2298	dtls(i, k) = dtls(i, k)/dtime
2299	dqls(i, k) = dqls(i, k)/dtime
2300	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
2301	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
2302	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
2303	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
2304	! c $ ,death_rate(i)*sigmaw(i)
2305	END IF
2306	END DO
2307	END DO
2308	!jyg<
2309	IF (iflag_wk_pop_dyn >= 1) THEN
2310	DO i = 1, klon
2311	IF (ok_qx_qw(i)) THEN
2312	d_sig_gen2(i) = d_sig_gen2(i)/dtime
2313	d_sig_death2(i) = d_sig_death2(i)/dtime
2314	d_sig_col2(i) = d_sig_col2(i)/dtime
2315	d_sig_spread2(i) = d_sig_spread2(i)/dtime
2316	d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
2317	d_sigmaw2(i) = d_sigmaw2(i)/dtime
2318	!
2319	d_dens_gen2(i) = d_dens_gen2(i)/dtime
2320	d_dens_death2(i) = d_dens_death2(i)/dtime
2321	d_dens_col2(i) = d_dens_col2(i)/dtime
2322	d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
2323	d_awdens2(i) = d_awdens2(i)/dtime
2324	d_wdens2(i) = d_wdens2(i)/dtime
2325	ENDIF
2326	ENDDO
2327	IF (iflag_wk_pop_dyn >= 2) THEN
2328	DO i = 1, klon
2329	IF (ok_qx_qw(i)) THEN
2330	d_adens_death2(i) = d_adens_death2(i)/dtime
2331	d_adens_icol2(i) = d_adens_icol2(i)/dtime
2332	d_adens_acol2(i) = d_adens_acol2(i)/dtime
2333	d_adens_bnd2(i) = d_adens_bnd2(i)/dtime
2334	ENDIF
2335	ENDDO
2336	IF (iflag_wk_pop_dyn == 3) THEN
2337	DO i = 1, klon
2338	IF (ok_qx_qw(i)) THEN
2339	d_asig_death2(i) = d_asig_death2(i)/dtime
2340	d_asig_iicol2(i) = d_asig_iicol2(i)/dtime
2341	d_asig_aicol2(i) = d_asig_aicol2(i)/dtime
2342	d_asig_spread2(i) = d_asig_spread2(i)/dtime
2343	d_asig_bnd2(i) = d_asig_bnd2(i)/dtime
2344	ENDIF
2345	ENDDO
2346	ENDIF ! (iflag_wk_pop_dyn == 3)
2347	ENDIF ! (iflag_wk_pop_dyn >= 2)
2348	ENDIF ! (iflag_wk_pop_dyn >= 1)
2349
2350	!>jyg
2351
2352	RETURN
2353	END SUBROUTINE wake
2354
2355	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, &
2356	d_deltaqw, sigmaw, d_sigmaw, alpha)
2357	! ------------------------------------------------------
2358	! Dtermination du coefficient alpha tel que les tendances
2359	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
2360	! a une humidite positive dans la zone (x) et dans la zone (w).
2361	! ------------------------------------------------------
2362	IMPLICIT NONE
2363
2364	! Input
2365	REAL qb(nlon, nl), d_qb(nlon, nl)
2366	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
2367	REAL sigmaw(nlon), d_sigmaw(nlon)
2368	LOGICAL wk_adv(nlon)
2369	INTEGER nl, nlon
2370	! Output
2371	REAL alpha(nlon)
2372	! Internal variables
2373	REAL zeta(nlon, nl)
2374	REAL alpha1(nlon)
2375	REAL x, a, b, c, discrim
2376	REAL epsilon_loc
2377	INTEGER i,k
2378
2379	DO k = 1, nl
2380	DO i = 1, nlon
2381	IF (wk_adv(i)) THEN
2382	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
2383	zeta(i, k) = 0.
2384	ELSE
2385	zeta(i, k) = 1.
2386	END IF
2387	END IF
2388	END DO
2389	DO i = 1, nlon
2390	IF (wk_adv(i)) THEN
2391	x = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qb(i, k) + &
2392	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i) &
2393	(deltaqw(i,k)+d_deltaqw(i,k))
2394	a = -d_sigmaw(i)*d_deltaqw(i, k)
2395	b = d_qb(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
2396	deltaqw(i, k)*d_sigmaw(i)
2397	c = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc
2398	discrim = bb - 4.a*c
2399	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
2400	IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap
2401	alpha1(i) = 1.
2402	ELSE
2403	IF (x>=0.) THEN
2404	alpha1(i) = 1.
2405	ELSE
2406	IF (a>0.) THEN
2407	alpha1(i) = 0.9min( (2.c)/(-b+sqrt(discrim)), &
2408	(-b+sqrt(discrim))/(2.*a) )
2409	ELSE IF (a==0.) THEN
2410	alpha1(i) = 0.9*(-c/b)
2411	ELSE
2412	! print*,'a,b,c discrim',a,b,c discrim
2413	alpha1(i) = 0.9max( (2.c)/(-b+sqrt(discrim)), &
2414	(-b+sqrt(discrim))/(2.*a))
2415	END IF
2416	END IF
2417	END IF
2418	alpha(i) = min(alpha(i), alpha1(i))
2419	END IF
2420	END DO
2421	END DO
2422
2423	RETURN
2424	END SUBROUTINE wake_vec_modulation
2425
2426
2427
2428	SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
2429	dth, hw_, rho, delta_t_min_in, &
2430	ktop, wk_adv, h_zzz, ptop1, ktop1)
2431
2432	USE lmdz_wake_ini , ONLY : wk_pupper
2433	USE lmdz_wake_ini , ONLY : RG
2434	USE lmdz_wake_ini , ONLY : hwmin
2435	USE lmdz_wake_ini , ONLY : iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t
2436	USE lmdz_wake_ini , ONLY : wk_int_delta_t_min
2437
2438	IMPLICIT NONE
2439
2440	INTEGER, INTENT(IN) :: klon,klev
2441	REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p
2442	REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: rho
2443	LOGICAL, DIMENSION (klon) , INTENT(IN) :: wk_adv
2444	REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: dth
2445	REAL, INTENT(IN) :: delta_t_min_in
2446
2447
2448	REAL, DIMENSION (klon) , INTENT(OUT) :: hw_
2449	REAL, DIMENSION (klon) , INTENT(OUT) :: ptop
2450	INTEGER, DIMENSION (klon) , INTENT(OUT) :: Ktop
2451	REAL, DIMENSION (klon) , INTENT(OUT) :: pupper
2452	INTEGER, DIMENSION (klon) , INTENT(OUT) :: kupper
2453	REAL, DIMENSION (klon) , INTENT(OUT) :: h_zzz !!
2454	REAL, DIMENSION (klon) , INTENT(OUT) :: Ptop1 !!
2455	INTEGER, DIMENSION (klon) , INTENT(OUT) :: ktop1 !!
2456
2457	INTEGER :: i,k
2458
2459	LOGICAL, DIMENSION (klon) :: wk_active
2460	REAL :: delta_t_min
2461	REAL, DIMENSION (klon) :: dthmin
2462	REAL, DIMENSION (klon) :: ptop_provis,ptop_new
2463	REAL, DIMENSION (klon) :: z, dz
2464	REAL, DIMENSION (klon) :: sum_dth
2465
2466	INTEGER, DIMENSION (klon) :: k_ptop_provis
2467	REAL, DIMENSION (klon) :: zk_ptop_provis
2468	REAL, DIMENSION (klon) :: omega !!
2469	REAL, DIMENSION (klon,klev+1) :: int_dth !!
2470	REAL, DIMENSION (klon,klev+1) :: dzz !!
2471	REAL, DIMENSION (klon,klev+1) :: zzz !!
2472	REAL, DIMENSION (klon) :: frac_int_dth !!
2473	REAL :: ddd!!
2474
2475
2476	INTEGER, SAVE :: ipas=0
2477	!$OMP THREADPRIVATE(ipas)
2478
2479
2480
2481	!INTEGER, SAVE :: compte=0
2482
2483	! LJYF : a priori z, dz sum_dth sont aussi des variables internes
2484	! Les eliminer apres verification convergence numerique
2485
2486	!compte=compte+1
2487	!print*,'compte=',compte
2488
2489	! Determine Ptop from buoyancy integral
2490	! ---------------------------------------
2491
2492	! - 1/ Pressure of the level where dth changes sign.
2493	!print*,'WAKE LJYF'
2494
2495
2496	if (iflag_wk_new_ptop==0) then
2497	delta_t_min=delta_t_min_in
2498	else
2499	delta_t_min=wk_delta_t_min
2500	endif
2501
2502	DO i = 1, klon
2503	ptop_provis(i) = ph(i, 1)
2504	k_ptop_provis(i) = 1
2505	END DO
2506
2507	DO k = 2, klev
2508	DO i = 1, klon
2509	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
2510	! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2511	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2512	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
2513	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
2514	k_ptop_provis(i) = k
2515	END IF
2516	END DO
2517	END DO
2518
2519
2520
2521	! - 2/ dth integral
2522
2523	DO i = 1, klon
2524	IF (wk_adv(i)) THEN !!! nrlmd
2525	sum_dth(i) = 0.
2526	dthmin(i) = -delta_t_min
2527	z(i) = 0.
2528	END IF
2529	END DO
2530
2531	DO k = 1, klev
2532	DO i = 1, klon
2533	IF (wk_adv(i)) THEN
2534	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
2535	IF (dz(i)>0) THEN
2536	z(i) = z(i) + dz(i)
2537	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
2538	dthmin(i) = amin1(dthmin(i), dth(i,k))
2539	END IF
2540	END IF
2541	END DO
2542	END DO
2543
2544	! - 3/ height of triangle with area= sum_dth and base = dthmin
2545
2546	DO i = 1, klon
2547	IF (wk_adv(i)) THEN
2548	hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
2549	hw_(i) = amax1(hwmin, hw_(i))
2550	END IF
2551	END DO
2552
2553	! - 4/ now, get Ptop
2554
2555	DO i = 1, klon
2556	IF (wk_adv(i)) THEN !!! nrlmd
2557	ktop(i) = 0
2558	z(i) = 0.
2559	END IF
2560	END DO
2561
2562	DO k = 1, klev
2563	DO i = 1, klon
2564	IF (wk_adv(i)) THEN
2565	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i))
2566	IF (dz(i)>0) THEN
2567	z(i) = z(i) + dz(i)
2568	ptop(i) = ph(i, k) - rho(i, k)RGdz(i)
2569	ktop(i) = k
2570	END IF
2571	END IF
2572	END DO
2573	END DO
2574
2575	! 4.5/Correct ktop and ptop
2576
2577	DO i = 1, klon
2578	ptop_new(i) = ptop(i)
2579	END DO
2580
2581	DO k = klev, 2, -1
2582	DO i = 1, klon
2583	! IM v3JYG; IF (k .GE. ktop(i)
2584	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
2585	! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2586	dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
2587	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
2588	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
2589	END IF
2590	END DO
2591	END DO
2592
2593
2594	DO i = 1, klon
2595	ptop(i) = ptop_new(i)
2596	END DO
2597
2598	DO k = klev, 1, -1
2599	DO i = 1, klon
2600	IF (wk_adv(i)) THEN !!! nrlmd
2601	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
2602	END IF
2603	END DO
2604	END DO
2605
2606	! IF (prt_level>=10) THEN
2607	! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
2608	! ENDIF
2609
2610	! -----------------------------------------------------------------------
2611	! nouveau calcul de hw et ptop
2612	! -----------------------------------------------------------------------
2613	!if (iflag_wk_new_ptop>0) then
2614	do i=1,klon
2615	ptop1(i)=ph(i,1)
2616	ktop1(i)=1
2617	h_zzz(i)=0.
2618	enddo
2619
2620	IF (iflag_wk_new_ptop/=0) THEN
2621
2622	int_dth(1:klon,1:klev+1)=0.
2623	DO i = 1, klon
2624	IF (wk_adv(i)) THEN
2625	int_dth(i,1) = 0.
2626	END IF
2627	END DO
2628
2629	if (abs(iflag_wk_new_ptop) == 1 ) then
2630	DO k = 2, klev+1
2631	Do i = 1, klon
2632	IF (wk_adv(i)) THEN
2633	if (k<=k_ptop_provis(i)) then
2634	ddd=dth(i,k-1)*(ph(i,k-1) - max(ptop_provis(i),ph(i,k)))
2635	!ddd=dth(i,k-1)*(ph(i,k-1) - ph(i,k))
2636	else
2637	ddd=0.
2638	endif
2639	int_dth(i,k) = int_dth(i,k-1) + ddd
2640	!ELSE
2641	! int_dth(i,k) = 0.
2642	END IF
2643	END DO
2644	END DO
2645	else
2646	k_ptop_provis(:)=klev+1
2647	dthmin(:)=dth(:,1)
2648	! calcul de l'int??grale de dT * dP jusqu'au dernier
2649	! niveau avec dT<0. (en s'assurant qu'on a bien un
2650	! dT negatif plus bas)
2651	DO k = 1, klev
2652	DO i = 1, klon
2653	dthmin(i)=min(dthmin(i),dth(i,k))
2654	ddd=dth(i,k)*(ph(i,k)-ph(i,k+1))
2655	if (dthmin(i)<0.) then
2656	if (k>=k_ptop_provis(i)) then
2657	ddd=0.
2658	else if (dth(i,k)>=0.) then
2659	ddd=0.
2660	k_ptop_provis(i)=k+1
2661	endif
2662	endif
2663	int_dth(i,k+1) = int_dth(i,k)+ ddd
2664	ENDDO
2665	ENDDO
2666
2667	DO i = 1, klon
2668	if ( k_ptop_provis(i)==klev+1 .or. .not. wk_adv(i)) then
2669	k_ptop_provis(i)=1
2670	endif
2671	ENDDO
2672	endif ! (abs(iflag_wk_new_ptop) == 1 )
2673	! print*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev)
2674	! print*, 'xxx, k_ptop_provis', k_ptop_provis(1)
2675
2676
2677
2678	! On se limite ?? des poches avec integrale dT * dp < -wk_int_delta_t_min
2679	do i=1,klon
2680	if (int_dth(i,k_ptop_provis(i)) > -wk_int_delta_t_min .or. k_ptop_provis(i)==1) then
2681	!if (1==0) then
2682	wk_active(i)=.false.
2683	ptop(i)=ph(i,1)
2684	ktop(i)=1
2685	hw_(i)=0.
2686	else
2687	wk_active(i)=wk_adv(i)
2688	endif
2689	enddo
2690
2691	DO i=1,klon
2692	IF (wk_active(i)) THEN
2693	frac_int_dth(i)=wk_frac_int_delta_t*int_dth(i,k_ptop_provis(i))
2694	ENDIF
2695	ENDDO
2696	DO k = 1,klev
2697	DO i =1, klon
2698	IF (wk_active(i)) THEN
2699	IF (int_dth(i,k)>=frac_int_dth(i)) THEN
2700	ktop1(i) = min(k, k_ptop_provis(i))
2701	!ktop1(i) = k
2702	!print*,ipas,'yyy ktop1= ',ktop1
2703	ENDIF
2704	ENDIF
2705	END DO
2706	END DO
2707	!print*, 'LAMINE'
2708
2709	DO i = 1, klon
2710	IF (wk_active(i)) THEN
2711	!print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1
2712	ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i))
2713	if (ddd==0.) then
2714	omega(i)=0.
2715	else
2716	omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd
2717	endif
2718	!! print*,'OMEGA ',omega(i)
2719	END IF
2720	END DO
2721
2722	!! print*, 'xxx'
2723	DO i = 1, klon
2724	IF (wk_active(i)) THEN
2725	! print*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', &
2726	! int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1)
2727	! print*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ', &
2728	!e omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1)
2729	ptop1(i) = min((1 - omega(i))ph(i,ktop1(i)) + omega(i)ph(i,ktop1(i)+1), ph(i,1))
2730	END IF
2731	END DO
2732
2733	DO i=1, klon
2734	IF (wk_active(i)) THEN
2735	zzz(i, 1) = 0
2736	END IF
2737	END DO
2738	DO k = 1, klev
2739	DO i = 1, klon
2740	IF (wk_active(i)) THEN
2741	dzz(i,k) = (ph(i,k) - ph(i,k+1))/(rho(i,k)*RG)
2742	zzz(i,k+1) = zzz(i,k) + dzz(i,k)
2743	END IF
2744	END DO
2745	END DO
2746
2747	DO i =1, klon
2748	IF (wk_active(i)) THEN
2749	h_zzz(i) = max((1- omega(i))zzz(i,ktop1(i)) + omega(i)zzz(i,ktop1(i)+1), hwmin)
2750	END IF
2751	END DO
2752
2753
2754	ENDIF ! (iflag_wk_new_ptop/=0)
2755
2756	!if (iflag_wk_new_ptop==2) then
2757	IF (iflag_wk_new_ptop>0) THEN
2758	do i=1,klon
2759	ptop(i)=ptop1(i)
2760	ktop(i)=ktop1(i)
2761	hw_(i)=h_zzz(i)
2762	enddo
2763
2764	!endif
2765	ENDIF
2766
2767	kupper = 0
2768
2769	IF (wk_pupper<1.) THEN
2770	! Choose an integration bound well above wake top
2771	! -----------------------------------------------------------------
2772
2773	! Pupper = 50000. ! melting level
2774	! Pupper = 60000.
2775	! Pupper = 80000. ! essais pour case_e
2776	DO i = 1, klon
2777	! pupper(i) = 0.6*ph(i, 1)
2778	pupper(i) = wk_pupper*ph(i, 1)
2779	pupper(i) = max(pupper(i), 45000.)
2780	! cc Pupper(i) = 60000.
2781	END DO
2782
2783	ELSE
2784	DO i=1, klon
2785	! pupper(i) = wk_pupperptop(i)+(1.-wk_pupper)ph(i, 1)
2786	! pupper(i) = min( wk_pupperptop(i)+(1.-wk_pupper)ph(i, 1) , ptop(i)-50.)
2787	pupper(i) = min( wk_pupperptop(i)+(1.-wk_pupper)ph(i, 1) , ptop(i)-5000.)
2788	END DO
2789	END IF
2790
2791	! -5/ Determination de kupper
2792
2793	DO k = klev, 1, -1
2794	DO i = 1, klon
2795	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
2796	END DO
2797	END DO
2798
2799	! On evite kupper = 1 et kupper = klev
2800	DO i = 1, klon
2801	kupper(i) = max(kupper(i), 2)
2802	kupper(i) = min(kupper(i), klev-1)
2803	END DO
2804	!---------- FIN nouveau calcul hw et ptop -------------------------------------
2805
2806	IF (iflag_wk_new_ptop==999) THEN
2807	DO i = 1, klon
2808	hw_(i)=0.
2809	ptop(i)=ph(i,1)
2810	Ktop(i)=1
2811	pupper(i)=ph(i,2)
2812	kupper(i)=2
2813	h_zzz(i)=0.
2814	Ptop1(i)=ph(i,1)
2815	ENDDO
2816	ENDIF
2817
2818	zk_ptop_provis=k_ptop_provis
2819
2820	RETURN
2821	END SUBROUTINE pkupper
2822
2823
2824	SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
2825	wdensmin, &
2826	dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
2827	d_awdens, d_wdens, d_sigmaw, &
2828	iflag_wk_act, wk_adv, cin, wape, &
2829	drdt, &
2830	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
2831	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
2832	d_wdens_targ, d_sigmaw_targ)
2833
2834
2835	USE lmdz_wake_ini , ONLY : wake_ini
2836	USE lmdz_wake_ini , ONLY : prt_level,RG
2837	USE lmdz_wake_ini , ONLY : stark, wdens_ref
2838	USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
2839	!! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
2840	USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
2841	USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
2842
2843	IMPLICIT NONE
2844
2845	INTEGER, INTENT(IN) :: klon,klev
2846	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
2847	REAL, INTENT(IN) :: dtime
2848	REAL, INTENT(IN) :: dtimesub
2849	REAL, INTENT(IN) :: wdensmin
2850	REAL, DIMENSION (klon), INTENT(IN) :: wgen
2851	REAL, DIMENSION (klon), INTENT(IN) :: wdens
2852	REAL, DIMENSION (klon), INTENT(IN) :: awdens
2853	REAL, DIMENSION (klon), INTENT(IN) :: sigmaw
2854	REAL, DIMENSION (klon), INTENT(IN) :: cstar
2855	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape
2856	REAL, DIMENSION (klon), INTENT(IN) :: f_shear
2857	INTEGER, INTENT(IN) :: iflag_wk_act
2858
2859
2860	!
2861
2862	! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields
2863	! computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
2864	REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk
2865	REAL, DIMENSION (klon), INTENT(OUT) :: gfl
2866	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_awdens, d_wdens
2867	REAL, DIMENSION (klon), INTENT(OUT) :: drdt
2868	! Some components of the tendencies of state variables
2869	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
2870	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_spread
2871	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
2872	REAL, INTENT(OUT) :: d_wdens_targ, d_sigmaw_targ
2873
2874
2875	REAL :: delta_t_min
2876	INTEGER :: i, k
2877	REAL :: wdens0
2878	! IM 080208
2879	LOGICAL, DIMENSION (klon) :: gwake
2880
2881	! Variables liees a la dynamique de population
2882	REAL, DIMENSION(klon) :: act
2883	REAL, DIMENSION(klon) :: tau_wk_inv
2884	REAL, DIMENSION(klon) :: wape1_act, wape2_act
2885	LOGICAL, DIMENSION (klon) :: kill_wake
2886	REAL :: drdt_pos
2887	REAL :: tau_wk_inv_min
2888
2889
2890
2891	IF (iflag_wk_act == 0) THEN
2892	act(:) = 0.
2893	ELSEIF (iflag_wk_act == 1) THEN
2894	act(:) = 1.
2895	ELSEIF (iflag_wk_act ==2) THEN
2896	DO i = 1, klon
2897	IF (wk_adv(i)) THEN
2898	wape1_act(i) = abs(cin(i))
2899	wape2_act(i) = 2.*wape1_act(i) + 1.
2900	act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
2901	ENDIF ! (wk_adv(i))
2902	ENDDO
2903	ENDIF ! (iflag_wk_act ==2)
2904
2905	DO i = 1, klon
2906	IF (wk_adv(i)) THEN
2907	rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
2908	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
2909	END IF
2910	END DO
2911
2912	DO i = 1, klon
2913	IF (wk_adv(i)) THEN
2914	!! tau_wk(i) = max(rad_wk(i)/(3.cstar(i))((cstar(i)/cstart)**1.5 - 1), 100.)
2915	tau_wk_inv(i) = max( (3.cstar(i))/(rad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
2916	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
2917	drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
2918	(1 + 2f_shear(i)(2.sigmaw(i)-aa0wdens(i)) - 2.*sigmaw(i))
2919	!! (1 - 2sigmaw(i)(1.-f_shear(i)))
2920	drdt_pos=max(drdt(i),0.)
2921
2922	!! d_wdens(i) = ( wgen(i)(1.+2.(sigmaw(i)-sigmad)) &
2923	!! - wdens(i)*tau_wk_inv_min &
2924	!! - 2.gfl(i)wdens(i)Cstar(i) )dtimesub
2925	!jyg+mlt<
2926	d_awdens(i) = ( wgen(i) - (1./tau_cv)(awdens(i) - act(i)wdens(i)) )*dtimesub
2927	d_dens_gen(i) = wgen(i)
2928	d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
2929	d_dens_col(i) = -2.wdens(i)gfl(i)*drdt_pos
2930	d_dens_gen(i) = d_dens_gen(i)*dtimesub
2931	d_dens_death(i) = d_dens_death(i)*dtimesub
2932	d_dens_col(i) = d_dens_col(i)*dtimesub
2933
2934	d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
2935	!! d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min - &
2936	!! 2.wdens(i)gfl(i)drdt_pos )dtimesub
2937	!>jyg+mlt
2938	!
2939	!jyg<
2940	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
2941	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
2942	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
2943	d_wdens(i) = d_wdens_targ
2944	!! d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
2945	!>jyg
2946
2947	!jyg+mlt<
2948	!! d_sigmaw(i) = ( (1.-2f_shear(i)sigmaw(i))(gfl(i)Cstar(i)+wgen(i)*sigmad/wdens(i)) &
2949	!! + 2.f_shear(i)wgen(i)sigmaw(i)*2/wdens(i) &
2950	!! - sigmaw(i)tau_wk_inv_min )dtimesub
2951	d_sig_gen(i) = wgen(i)*aa0
2952	d_sig_death(i) = - sigmaw(i)(1.-awdens(i)/wdens(i))tau_wk_inv_min
2953	!!
2954
2955	d_sig_col(i) = - 2f_shear(i)sigmaw(i)gfl(i)drdt_pos
2956	d_sig_col(i) = - 2f_shear(i)(2.sigmaw(i)-wdens(i)aa0)gfl(i)drdt_pos
2957	d_sig_spread(i) = gfl(i)*cstar(i)
2958	d_sig_gen(i) = d_sig_gen(i)*dtimesub
2959	d_sig_death(i) = d_sig_death(i)*dtimesub
2960	d_sig_col(i) = d_sig_col(i)*dtimesub
2961	d_sig_spread(i) = d_sig_spread(i)*dtimesub
2962	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
2963	!>jyg+mlt
2964	!
2965	!jyg<
2966	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
2967	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
2968	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
2969	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
2970	d_sigmaw(i) = d_sigmaw_targ
2971	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
2972	!>jyg
2973	ENDIF
2974	ENDDO
2975
2976	IF (prt_level >= 10) THEN
2977	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
2978	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
2979	print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
2980	wdens(1), awdens(1), act(1), d_awdens(1)
2981	print ,'wake, wgen, -(wdens-awdens)tau_wk_inv, -2.wdensgfl*drdt_pos, d_wdens ', &
2982	wgen(1), -(wdens(1)-awdens(1))tau_wk_inv(1), -2.wdens(1)gfl(1)drdt_pos, d_wdens(1)
2983	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
2984	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
2985	ENDIF
2986
2987	RETURN
2988	END SUBROUTINE wake_popdyn_1
2989
2990	SUBROUTINE wake_popdyn_2 ( klon, klev, wk_adv, dtimesub, wgen, &
2991	wdensmin, &
2992	sigmaw, wdens, awdens, & !! states variables
2993	gfl, cstar, cin, wape, rad_wk, &
2994	d_sigmaw, d_wdens, d_awdens, & !! tendences
2995	cont_fact, &
2996	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
2997	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
2998	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
2999
3000
3001
3002	USE lmdz_wake_ini , ONLY : wake_ini
3003	USE lmdz_wake_ini , ONLY : prt_level,RG
3004	USE lmdz_wake_ini , ONLY : stark, wdens_ref
3005	USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
3006	!! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
3007	USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
3008	USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
3009
3010	IMPLICIT NONE
3011
3012	INTEGER, INTENT(IN) :: klon,klev
3013	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
3014	REAL, INTENT(IN) :: dtimesub
3015	REAL, INTENT(IN) :: wdensmin
3016	REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes
3017	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes
3018	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area
3019	REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area
3020	REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes
3021	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre
3022
3023	!
3024	REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = wake radius
3025	REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front lenght per unit area
3026	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_wdens, d_awdens
3027	REAL, DIMENSION (klon), INTENT(OUT) :: cont_fact !! RM facteur de contact = 2 pi * rad * C*
3028	! Some components of the tendencies of state variables
3029	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
3030	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
3031	REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd
3032
3033
3034	!! internal variables
3035
3036	INTEGER :: i, k
3037	REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes
3038	REAL :: tau_wk_inv_min
3039	REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes
3040	REAL :: d_wdens_targ, d_sigmaw_targ
3041
3042
3043	!! Equations
3044	!! dD/dt = B - (D-A)/tau - f D^2
3045	!! dA/dt = B - A/tau_prime + f (D-A)^2 - f A^2
3046	!! dsigma/dt = B a0 - sigma/D (D-A)/tau + Lg C* - f (D-A)^2 (sigma/D-a0)
3047	!!
3048	!! f = 2 (B (a0-sigma/D) + Lg C*) / (2 (D-A)^2 (2 sigma/D-a0) + D (1-2 sigma))
3049
3050
3051	DO i = 1, klon
3052	IF (wk_adv(i)) THEN
3053	rad_wk(i) = max( sqrt(sigmaw(i)/(3.14*wdens(i))) , rzero)
3054	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
3055	END IF
3056	END DO
3057
3058
3059	DO i = 1, klon
3060	IF (wk_adv(i)) THEN
3061	!! tau_wk(i) = max(rad_wk(i)/(3.cstar(i))((cstar(i)/cstart)**1.5 - 1), 100.)
3062	tau_wk_inv(i) = max( (3.cstar(i))/(rad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
3063	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
3064	tau_prime(i) = tau_cv
3065	!! cont_fact(i) = 2.(wgen(i)(aa0-sigmaw(i)/wdens(i)) + gfl(i)*cstar(i)) / &
3066	!! (2.(wdens(i)-awdens(i))2(2.sigmaw(i)/wdens(i) - aa0) + wdens(i)(1.-2.*sigmaw(i)))
3067	!! cont_fact(i) = 2.3.14rad_wk(i)*cstar(i) ! bug
3068	!! cont_fact(i) = 4.3.14rad_wk(i)*cstar(i)
3069	cont_fact(i) = 2.gfl(i)cstar(i)/wdens(i)
3070
3071	d_sig_gen(i) = wgen(i)*aa0
3072	d_sig_death(i) = - sigmaw(i)(1.-awdens(i)/wdens(i))tau_wk_inv_min
3073	d_sig_col(i) = - cont_fact(i)(wdens(i)-awdens(i))2(2.*sigmaw(i)/wdens(i)-aa0)
3074	d_sig_spread(i) = gfl(i)*cstar(i)
3075	!
3076	d_sig_gen(i) = d_sig_gen(i)*dtimesub
3077	d_sig_death(i) = d_sig_death(i)*dtimesub
3078	d_sig_col(i) = d_sig_col(i)*dtimesub
3079	d_sig_spread(i) = d_sig_spread(i)*dtimesub
3080	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
3081
3082
3083	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
3084	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
3085	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
3086	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
3087	d_sigmaw(i) = d_sigmaw_targ
3088	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
3089
3090
3091	d_dens_gen(i) = wgen(i)
3092	d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min
3093	d_dens_col(i) = - cont_fact(i)wdens(i)*2
3094	!
3095	d_dens_gen(i) = d_dens_gen(i)*dtimesub
3096	d_dens_death(i) = d_dens_death(i)*dtimesub
3097	d_dens_col(i) = d_dens_col(i)*dtimesub
3098	d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
3099
3100
3101	d_adens_death(i) = -awdens(i)/tau_prime(i)
3102	d_adens_icol(i) = cont_fact(i)(wdens(i)-awdens(i))*2
3103	d_adens_acol(i) = - cont_fact(i)awdens(i)*2
3104	!
3105	d_adens_death(i) = d_adens_death(i)*dtimesub
3106	d_adens_icol(i) = d_adens_icol(i)*dtimesub
3107	d_adens_acol(i) = d_adens_acol(i)*dtimesub
3108	d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)
3109
3110	!!
3111	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
3112	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
3113	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
3114	d_wdens(i) = d_wdens_targ
3115
3116	d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
3117	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
3118	d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
3119	d_awdens(i) = d_wdens_targ
3120
3121
3122
3123	ENDIF
3124	ENDDO
3125
3126	IF (prt_level >= 10) THEN
3127	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1) ', &
3128	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), cont_fact(1)
3129	print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
3130	wdens(1), awdens(1), d_awdens(1)
3131	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
3132	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
3133	ENDIF
3134	sigmaw=sigmaw+d_sigmaw
3135	wdens=wdens+d_wdens
3136	awdens=awdens+d_awdens
3137
3138	RETURN
3139	END SUBROUTINE wake_popdyn_2
3140
3141	SUBROUTINE wake_popdyn_3 ( klon, klev, phys_sub, wk_adv, dtimesub, wgen, &
3142	wdensmin, &
3143	sigmaw, asigmaw, wdens, awdens, & !! state variables
3144	gfl, agfl, cstar, cin, wape, &
3145	rad_wk, arad_wk, irad_wk, &
3146	d_sigmaw, d_asigmaw, d_wdens, d_awdens, & !! tendencies
3147	d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
3148	d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd, &
3149	d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
3150	d_adens_death, d_adens_icol, d_adens_acol, d_adens_bnd )
3151
3152
3153
3154	USE lmdz_wake_ini , ONLY : wake_ini
3155	USE lmdz_wake_ini , ONLY : prt_level,RG
3156	USE lmdz_wake_ini , ONLY : stark, wdens_ref
3157	USE lmdz_wake_ini , ONLY : tau_cv, rzero, aa0
3158	!! USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn, wdensmin
3159	USE lmdz_wake_ini , ONLY : iflag_wk_pop_dyn
3160	USE lmdz_wake_ini , ONLY : sigmad, cstart, sigmaw_max
3161	USE lmdz_wake_ini , ONLY : smallestreal
3162
3163	IMPLICIT NONE
3164
3165	INTEGER, INTENT(IN) :: klon,klev
3166	LOGICAL, INTENT(IN) :: phys_sub
3167	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
3168	REAL, INTENT(IN) :: dtimesub
3169	REAL, INTENT(IN) :: wdensmin
3170	REAL, DIMENSION (klon), INTENT(IN) :: wgen !! B = birth rate of wakes
3171	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw !! sigma = fractional area of wakes
3172	REAL, DIMENSION (klon), INTENT(INOUT) :: asigmaw !! sigma = fractional area of active wakes
3173	REAL, DIMENSION (klon), INTENT(INOUT) :: wdens !! D = number of wakes per unit area
3174	REAL, DIMENSION (klon), INTENT(INOUT) :: awdens !! A = number of active wakes per unit area
3175	REAL, DIMENSION (klon), INTENT(IN) :: cstar !! C* = spreading velocity of wakes
3176	REAL, DIMENSION (klon), INTENT(IN) :: cin, wape ! RM : A Faire disparaitre
3177
3178	!
3179	REAL, DIMENSION (klon), INTENT(OUT) :: rad_wk !! r = mean wake radius
3180	REAL, DIMENSION (klon), INTENT(OUT) :: arad_wk !! r_A = wake radius of active wakes
3181	REAL, DIMENSION (klon), INTENT(OUT) :: irad_wk !! r_I = wake radius of inactive wakes
3182	REAL, DIMENSION (klon), INTENT(OUT) :: gfl !! Lg = gust front length per unit area
3183	REAL, DIMENSION (klon), INTENT(OUT) :: agfl !! LgA = gust front length of active wakes
3184	!! per unit area
3185	REAL, DIMENSION (klon), INTENT(OUT) :: d_sigmaw, d_asigmaw, d_wdens, d_awdens
3186	! Some components of the tendencies of state variables
3187	REAL, DIMENSION (klon), INTENT(OUT) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
3188	REAL, DIMENSION (klon), INTENT(OUT) :: d_asig_death, d_asig_aicol, d_asig_iicol, d_asig_spread, d_asig_bnd
3189	REAL, DIMENSION (klon), INTENT(OUT) :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
3190	REAL, DIMENSION (klon), INTENT(OUT) :: d_adens_death, d_adens_acol, d_adens_icol, d_adens_bnd
3191
3192
3193	!! internal variables
3194
3195	INTEGER :: i, k
3196	REAL, DIMENSION (klon) :: iwdens, isigmaw !! inactive wake density and fractional area
3197	!! REAL, DIMENSION (klon) :: d_arad, d_irad
3198	REAL, DIMENSION (klon) :: igfl !! LgI = gust front length of inactive wakes
3199	!! per unit area
3200	REAL, DIMENSION (klon) :: s_wk !! mean area of individual wakes
3201	REAL, DIMENSION (klon) :: as_wk !! mean area of individual active wakes
3202	REAL, DIMENSION (klon) :: is_wk !! mean area of individual inactive wakes
3203	REAL, DIMENSION (klon) :: tau_wk_inv !! tau = life time of wakes
3204	REAL :: tau_wk_inv_min
3205	REAL, DIMENSION (klon) :: tau_prime !! tau_prime = life time of actives wakes
3206	REAL :: d_wdens_targ, d_sigmaw_targ
3207
3208
3209	!! Equations
3210	!! ---------
3211	!! Gust fronts:
3212	!! Lg_A = 2 pi r_A A
3213	!! Lg_I = 2 pi r_I I
3214	!! Lg = 2 pi r D
3215	!!
3216	!! Areas:
3217	!! s = pi r^2
3218	!! s_A = pi r_A^2
3219	!! s_I = pi r_I^2
3220	!!
3221	!! Life expectancy:
3222	!! tau_I = 3 C* ((C/Ct)^3/2 - 1) / r_I
3223	!!
3224	!! Time deratives:
3225	!! dD/dt = B - (D-A)/tau_I - 2 Lg C* D
3226	!! dA/dt = B - A/tau_A + 2 Lg_I C* (D-A) - 2 Lg_A C* A
3227	!! dsigma/dt = B a0 - sigma_I/tau_I + Lg C* - 2 Lg_I C* (D-A) (2 s_I - a0)
3228	!! dsigma_A/dt = B a0 - sigma_A/tau_A + Lg_A C* + (Lg_A I + Lg_I A) C* s_I + 2 Lg_I C* I a0
3229	!!
3230
3231	DO i = 1, klon
3232	IF (wk_adv(i)) THEN
3233	iwdens(i) = wdens(i) - awdens(i)
3234	isigmaw(i) = sigmaw(i) - asigmaw(i)
3235	!
3236	arad_wk(i) = max( sqrt(asigmaw(i)/(3.14*awdens(i))) , rzero)
3237	irad_wk(i) = max( sqrt((sigmaw(i)-asigmaw(i))/ &
3238	(3.14*max(smallestreal,(wdens(i)-awdens(i))))), rzero)
3239	rad_wk(i) = (awdens(i)arad_wk(i)+(wdens(i)-awdens(i))irad_wk(i))/wdens(i)
3240	!
3241	s_wk(i) = 3.14rad_wk(i)*2
3242	as_wk(i) = 3.14arad_wk(i)*2
3243	is_wk(i) = 3.14irad_wk(i)*2
3244	!
3245	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
3246	agfl(i) = 2.sqrt(3.14awdens(i)*asigmaw(i))
3247	igfl(i) = gfl(i) - agfl(i)
3248	ENDIF
3249	ENDDO
3250
3251
3252	DO i = 1, klon
3253	IF (wk_adv(i)) THEN
3254	tau_wk_inv(i) = max( (3.cstar(i))/(irad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
3255	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
3256	tau_prime(i) = tau_cv
3257
3258	d_sig_gen(i) = wgen(i)*aa0
3259	d_sig_death(i) = - isigmaw(i)*tau_wk_inv_min
3260	d_sig_col(i) = - 2.igfl(i)cstar(i)iwdens(i)(2.*is_wk(i)-aa0)
3261	d_sig_spread(i) = gfl(i)*cstar(i)
3262	!
3263	d_sig_gen(i) = d_sig_gen(i)*dtimesub
3264	d_sig_death(i) = d_sig_death(i)*dtimesub
3265	d_sig_col(i) = d_sig_col(i)*dtimesub
3266	d_sig_spread(i) = d_sig_spread(i)*dtimesub
3267	d_sigmaw(i) = d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
3268	IF (CPPKEY_IOPHYS_WK) THEN
3269	IF (phys_sub) call iophys_ecrit('d_sigmaw0',1,'d_sigmaw0','',d_sigmaw)
3270	END IF
3271
3272
3273	d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
3274	!! d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
3275	!! d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
3276	d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
3277	d_sigmaw(i) = d_sigmaw_targ
3278	!! d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
3279	IF (CPPKEY_IOPHYS_WK) THEN
3280	IF (phys_sub) THEN
3281	call iophys_ecrit('tauwk_inv',1,'tau_wk_inv_min','',tau_wk_inv_min)
3282	call iophys_ecrit('d_sigmaw',1,'d_sigmaw','',d_sigmaw)
3283	call iophys_ecrit('d_sig_gen',1,'d_sig_gen','',d_sig_gen)
3284	call iophys_ecrit('d_sig_death',1,'d_sig_death','',d_sig_death)
3285	call iophys_ecrit('d_sig_col',1,'d_sig_col','',d_sig_col)
3286	call iophys_ecrit('d_sig_spread',1,'d_sig_spread','',d_sig_spread)
3287	call iophys_ecrit('d_sig_bnd',1,'d_sig_bnd','',d_sig_bnd)
3288	ENDIF
3289	END IF
3290	d_asig_death(i) = - asigmaw(i)/tau_prime(i)
3291	d_asig_aicol(i) = (agfl(i)iwdens(i) + igfl(i)awdens(i))cstar(i)is_wk(i)
3292	d_asig_iicol(i) = 2.igfl(i)cstar(i)iwdens(i)aa0
3293	d_asig_spread(i) = agfl(i)*cstar(i)
3294	!
3295	d_asig_death(i) = d_asig_death(i)*dtimesub
3296	d_asig_aicol(i) = d_asig_aicol(i)*dtimesub
3297	d_asig_iicol(i) = d_asig_iicol(i)*dtimesub
3298	d_asig_spread(i) = d_asig_spread(i)*dtimesub
3299	d_asigmaw(i) = d_sig_gen(i) + d_asig_death(i) + d_asig_aicol(i) + d_asig_iicol(i) + d_asig_spread(i)
3300	IF (CPPKEY_IOPHYS_WK) THEN
3301	IF (phys_sub) call iophys_ecrit('d_asigmaw0',1,'d_asigmaw0','',d_asigmaw)
3302	END IF
3303
3304	d_sigmaw_targ = min(max(d_asigmaw(i),-asigmaw(i)), sigmaw(i)-asigmaw(i))
3305	!! d_dens_bnd(i) = d_dens_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
3306	d_asig_bnd(i) = d_sigmaw_targ - d_asigmaw(i)
3307	d_asigmaw(i) = d_sigmaw_targ
3308	IF (CPPKEY_IOPHYS_WK) THEN
3309	IF (phys_sub) THEN
3310	call iophys_ecrit('d_asigmaw',1,'d_asigmaw','',d_asigmaw)
3311	call iophys_ecrit('d_asig_death',1,'d_asig_death','',d_asig_death)
3312	call iophys_ecrit('d_asig_aicol',1,'d_asig_aicol','',d_asig_aicol)
3313	call iophys_ecrit('d_asig_iicol',1,'d_asig_iicol','',d_asig_iicol)
3314	call iophys_ecrit('d_asig_spread',1,'d_asig_spread','',d_asig_spread)
3315	call iophys_ecrit('d_asig_bnd',1,'d_asig_bnd','',d_asig_bnd)
3316	ENDIF
3317	END IF
3318	d_dens_gen(i) = wgen(i)
3319	d_dens_death(i) = - iwdens(i)*tau_wk_inv_min
3320	d_dens_col(i) = - 2.gfl(i)cstar(i)*wdens(i)
3321	!
3322	d_dens_gen(i) = d_dens_gen(i)*dtimesub
3323	d_dens_death(i) = d_dens_death(i)*dtimesub
3324	d_dens_col(i) = d_dens_col(i)*dtimesub
3325	d_wdens(i) = d_dens_gen(i) + d_dens_death(i) + d_dens_col(i)
3326	!!
3327	d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
3328	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
3329	d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
3330	d_wdens(i) = d_wdens_targ
3331	IF (CPPKEY_IOPHYS_WK) THEN
3332	IF (phys_sub) THEN
3333	call iophys_ecrit('d_wdens',1,'d_wdens','',d_wdens)
3334	call iophys_ecrit('d_dens_gen',1,'d_dens_gen','',d_dens_gen)
3335	call iophys_ecrit('d_dens_death',1,'d_dens_death','',d_dens_death)
3336	call iophys_ecrit('d_dens_col',1,'d_dens_col','',d_dens_col)
3337	ENDIF
3338	END IF
3339
3340	d_adens_death(i) = -awdens(i)/tau_prime(i)
3341	d_adens_icol(i) = 2.igfl(i)cstar(i)*iwdens(i)
3342	d_adens_acol(i) = - 2.agfl(i)cstar(i)*awdens(i)
3343	!
3344	d_adens_death(i) = d_adens_death(i)*dtimesub
3345	d_adens_icol(i) = d_adens_icol(i)*dtimesub
3346	d_adens_acol(i) = d_adens_acol(i)*dtimesub
3347	d_awdens(i) = d_dens_gen(i) + d_adens_death(i) + d_adens_icol(i) + d_adens_acol(i)
3348	IF (CPPKEY_IOPHYS_WK) THEN
3349	IF (phys_sub) THEN
3350	call iophys_ecrit('d_awdens',1,'d_awdens','',d_awdens)
3351	call iophys_ecrit('d_adens_death',1,'d_adens_death','',d_adens_death)
3352	call iophys_ecrit('d_adens_icol',1,'d_adens_icol','',d_adens_icol)
3353	call iophys_ecrit('d_adens_acol',1,'d_adens_acol','',d_adens_acol)
3354	ENDIF
3355	END IF
3356	d_wdens_targ = min(max(d_awdens(i),-awdens(i)), wdens(i)-awdens(i))
3357	!! d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
3358	d_adens_bnd(i) = d_wdens_targ - d_awdens(i)
3359	d_awdens(i) = d_wdens_targ
3360
3361	!! d_irad(i) = (d_sigmaw(i)-d_asigmaw(i)-isigmaw(i)*(d_wdens(i)-awdens(i))/iwdens(i)) / &
3362	!! max(smallestreal,(2.3.14iwdens(i)*irad_wk(i)))
3363	!! d_arad(i) = (d_asigmaw(i)-asigmaw(i)*d_awdens(i)/awdens(i)) / &
3364	!! max(smallestreal,(2.3.14awdens(i)*arad_wk(i)))
3365	!! d_irad(i) = d_irad(i)*dtimesub
3366	!! d_arad(i) = d_arad(i)*dtimesub
3367	!! call iophys_ecrit('d_irad',1,'d_irad','m',d_irad)
3368	!! call iophys_ecrit('d_airad',1,'d_arad','m',d_arad)
3369	!!
3370	ENDIF
3371	ENDDO
3372
3373	IF (prt_level >= 10) THEN
3374	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1) ', &
3375	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), gfl(1)
3376	print *,'wake, wdens(1), awdens(1), d_awdens(1) ', &
3377	wdens(1), awdens(1), d_awdens(1)
3378	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
3379	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
3380	ENDIF
3381	sigmaw=sigmaw+d_sigmaw
3382	asigmaw=asigmaw+d_asigmaw
3383	wdens=wdens+d_wdens
3384	awdens=awdens+d_awdens
3385
3386	RETURN
3387	END SUBROUTINE wake_popdyn_3
3388
3389	END MODULE lmdz_wake

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

Download in other formats:

Original Format