Context Navigation

wake.F90 @ 4244

Last change on this file since 4244 was 4231, checked in by fhourdin, 2 years ago
Removing prints in wake.F90
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: 85.9 KB

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

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

Download in other formats:

Original Format