Context Navigation

lmdz_wake.F90 @ 5212

Last change on this file since 5212 was 5203, checked in by Laurent Fairhead, 7 weeks ago
Merge with trunk revision 5202 before reintegration to trunk
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.3 KB

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

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

Download in other formats:

Original Format