Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lmdz_wake.F90 @ 5213

Last change on this file since 5213 was 4727, checked in by idelkadi, 13 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
File size: 82.1 KB

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

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

Download in other formats:

Original Format