Context Navigation

lmdz_wake.f90 @ 5506

Last change on this file since 5506 was 5499, checked in by yann meurdesoif, 3 days ago
Adapt calwake and wake for automatic GPU port. YM
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.4 KB

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

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

Download in other formats:

Original Format