Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lmdz_wake.F90 @ 5451

Last change on this file since 5451 was 4727, checked in by idelkadi, 15 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
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.8 KB

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

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

Download in other formats:

Original Format