Context Navigation

lmdz_wake.F90 @ 5418

Last change on this file since 5418 was 5285, checked in by abarral, 7 weeks ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 81.9 KB

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

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

Download in other formats:

Original Format