Context Navigation

wake.F90 @ 3627

Last change on this file since 3627 was 3605, checked in by lguez, 5 years ago
Merge revisions 3427:3600 of trunk into branch `Ocean_skin`
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 71.7 KB

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

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

Download in other formats:

Original Format