Context Navigation

lmdz_wake.F90 @ 5006

Last change on this file since 5006 was 4908, checked in by jyg, 9 months ago
New wake height formulation in the wake parametrization
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.1 KB

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

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

Download in other formats:

Original Format