Context Navigation

lmdz_wake.F90 @ 5218

Last change on this file since 5218 was 5193, checked in by abarral, 3 months ago
Merge r5181, r5188 from trunk Lint lmdz_wake.F90, lmdz_wake_ini.F90
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 120.7 KB

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

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

Download in other formats:

Original Format