Context Navigation

wake.F90 @ 5454

Last change on this file since 5454 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
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: 88.8 KB

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

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

Download in other formats:

Original Format