Context Navigation

lmdz_wake.F90 @ 4670

Last change on this file since 4670 was 4588, checked in by fhourdin, 13 months ago
Nouveaux noms et encapsulation module pour wake
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: 89.2 KB

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

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

Download in other formats:

Original Format