Context Navigation

wake.F90 @ 3999

Last change on this file since 3999 was 3648, checked in by jghattas, 5 years ago
Correction d'erreur introduit en commit [3641] pour faire tourner en mode debug.
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: 72.2 KB

Rev	Line
[1992]	1
[1403]	2	! $Id: wake.F90 3648 2020-03-16 15:36:59Z evignon $
[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
[3455]	198	REAL, SAVE :: wdensmin
	199	!$OMP THREADPRIVATE(wdensmin)
[3208]	200	REAL, SAVE :: sigmad, hwmin, wapecut, cstart
	201	!$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
[3455]	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./
[3641]	319	!! DATA wapecut, sigmad, hwmin/1., .02, 10./
	320	DATA sigmad, hwmin/.02, 10./
[3252]	321	!! DATA wdensmin/1.e-12/
	322	DATA wdensmin/1.e-14/
[1992]	323	! cc nrlmd
	324	DATA sigmaw_max/0.4/
	325	DATA dens_rate/0.1/
	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
[3648]	362	! Get wapecut from parameter file
	363	wapecut = 1.
	364	CALL getin_p('wapecut', wapecut)
	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
[3641]	381	!
	382	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	383	!!!!!!!!! Population dynamics parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!
	384	!------------------------------------------------------------------------
	385
[3208]	386	iflag_wk_pop_dyn = 0
	387	CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed
	388	! and wdens prognostic
	389	iflag_wk_act = 0
	390	CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
	391	! 1: act(:)=1.
	392	! 2: act(:)=f(Wape)
[3648]	393
[3641]	394	rzero = 5000.
	395	CALL getin_p('rzero_wk', rzero)
	396	aa0 = 3.14rzerorzero
	397	!
	398	tau_cv = 4000.
	399	CALL getin_p('tau_cv', tau_cv)
	400
	401	!------------------------------------------------------------------------
	402
[2474]	403	coefgw=4.
	404	CALL getin_p('coefgw',coefgw)
	405
	406	WRITE(,) 'stark=', stark
	407	WRITE(,) 'alpk=', alpk
[2761]	408	!jyg<
	409	!! WRITE(,) 'wdens_ref=', wdens_ref
	410	WRITE(,) 'wdens_ref_o=', wdens_ref(1)
	411	WRITE(,) 'wdens_ref_l=', wdens_ref(2)
	412	!>jyg
[3208]	413	WRITE(,) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
	414	WRITE(,) 'iflag_wk_act',iflag_wk_act
[2474]	415	WRITE(,) 'coefgw=', coefgw
	416
[2757]	417	flag_wk_check_trgl=.false.
	418	CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
	419	WRITE(,) 'flag_wk_check_trgl=', flag_wk_check_trgl
[2922]	420	WRITE(,) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
	421	iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
	422	CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
	423	WRITE(,) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
[2757]	424
[2078]	425	first=.false.
	426	endif
	427
[3208]	428	IF (iflag_wk_pop_dyn == 0) THEN
[1992]	429	! Initialisation de toutes des densites a wdens_ref.
	430	! Les densites peuvent evoluer si les poches debordent
	431	! (voir au tout debut de la boucle sur les substeps)
[3208]	432	!jyg<
	433	!! wdens(:) = wdens_ref
	434	DO i = 1,klon
	435	wdens(i) = wdens_ref(znatsurf(i)+1)
	436	ENDDO
	437	!>jyg
	438	ENDIF ! (iflag_wk_pop_dyn == 0)
[974]	439
[1992]	440	! print*,'stark',stark
	441	! print*,'alpk',alpk
	442	! print*,'wdens',wdens
	443	! print*,'coefgw',coefgw
	444	! cc
	445	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
	446	! -------------------------------------------------------------------------
[974]	447
[1992]	448	delta_t_min = 0.2
[974]	449
[2671]	450	! 1. - Save initial values, initialize tendencies, initialize output fields
	451	! ------------------------------------------------------------------------
[974]	452
[2671]	453	!jyg<
	454	!! DO k = 1, klev
	455	!! DO i = 1, klon
	456	!! ppi(i, k) = pi(i, k)
	457	!! deltatw0(i, k) = deltatw(i, k)
	458	!! deltaqw0(i, k) = deltaqw(i, k)
	459	!! te(i, k) = te0(i, k)
	460	!! qe(i, k) = qe0(i, k)
	461	!! dtls(i, k) = 0.
	462	!! dqls(i, k) = 0.
	463	!! d_deltat_gw(i, k) = 0.
	464	!! d_te(i, k) = 0.
	465	!! d_qe(i, k) = 0.
	466	!! d_deltatw(i, k) = 0.
	467	!! d_deltaqw(i, k) = 0.
	468	!! ! IM 060508 beg
	469	!! d_deltatw2(i, k) = 0.
	470	!! d_deltaqw2(i, k) = 0.
	471	!! ! IM 060508 end
	472	!! END DO
	473	!! END DO
	474	ppi(:,:) = pi(:,:)
	475	deltatw0(:,:) = deltatw(:,:)
	476	deltaqw0(:,:) = deltaqw(:,:)
	477	te(:,:) = te0(:,:)
	478	qe(:,:) = qe0(:,:)
	479	dtls(:,:) = 0.
	480	dqls(:,:) = 0.
	481	d_deltat_gw(:,:) = 0.
	482	d_te(:,:) = 0.
	483	d_qe(:,:) = 0.
	484	d_deltatw(:,:) = 0.
	485	d_deltaqw(:,:) = 0.
	486	d_deltatw2(:,:) = 0.
	487	d_deltaqw2(:,:) = 0.
[3208]	488
	489	IF (iflag_wk_act == 0) THEN
	490	act(:) = 0.
	491	ELSEIF (iflag_wk_act == 1) THEN
	492	act(:) = 1.
	493	ENDIF
	494
[2671]	495	!! DO i = 1, klon
	496	!! sigmaw_in(i) = sigmaw(i)
	497	!! END DO
	498	sigmaw_in(:) = sigmaw(:)
	499	!>jyg
	500
[1992]	501	! sigmaw1=sigmaw
	502	! IF (sigd_con.GT.sigmaw1) THEN
	503	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
	504	! ENDIF
[3208]	505	IF (iflag_wk_pop_dyn >=1) THEN
	506	DO i = 1, klon
	507	wdens_targ = max(wdens(i),wdensmin)
	508	d_wdens2(i) = wdens_targ - wdens(i)
	509	wdens(i) = wdens_targ
	510	END DO
	511	ELSE
	512	DO i = 1, klon
	513	d_awdens2(i) = 0.
	514	d_wdens2(i) = 0.
	515	END DO
	516	ENDIF ! (iflag_wk_pop_dyn >=1)
	517	!
[1992]	518	DO i = 1, klon
	519	! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
[2635]	520	!jyg<
	521	!! sigmaw(i) = amax1(sigmaw(i), sigmad)
	522	!! sigmaw(i) = amin1(sigmaw(i), 0.99)
	523	sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
	524	d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
	525	sigmaw(i) = sigmaw_targ
	526	!>jyg
[1992]	527	END DO
[3208]	528
[2671]	529	!
[3208]	530	IF (iflag_wk_pop_dyn >= 1) THEN
	531	awdens_in(:) = awdens(:)
	532	wdens_in(:) = wdens(:)
	533	!! wdens(:) = wdens(:) + wgen(:)*dtime
	534	!! d_wdens2(:) = wgen(:)*dtime
	535	!! ELSE
	536	ENDIF ! (iflag_wk_pop_dyn >= 1)
	537
	538	wape(:) = 0.
	539	wape2(:) = 0.
	540	d_sigmaw(:) = 0.
	541	ktopw(:) = 0
	542	!
[2671]	543	!<jyg
	544	dth(:,:) = 0.
	545	tu(:,:) = 0.
	546	qu(:,:) = 0.
	547	dtke(:,:) = 0.
	548	dqke(:,:) = 0.
	549	spread(:,:) = 0.
	550	omgbdth(:,:) = 0.
	551	omg(:,:) = 0.
	552	dp_omgb(:,:) = 0.
	553	dp_deltomg(:,:) = 0.
	554	hw(:) = 0.
	555	wape(:) = 0.
	556	fip(:) = 0.
	557	gfl(:) = 0.
	558	cstar(:) = 0.
	559	ktopw(:) = 0
	560	!
	561	! Vertical advection local variables
	562	omgbw(:,:) = 0.
	563	omgtop(:) = 0
	564	dp_omgbw(:,:) = 0.
	565	omgbdq(:,:) = 0.
	566	!>jyg
	567	!
	568	IF (prt_level>=10) THEN
	569	PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
	570	PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
	571	PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
	572	PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
	573	PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
	574	PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
	575	PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
	576	PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
	577	PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
	578	ENDIF
[974]	579
[1992]	580	! 2. - Prognostic part
	581	! --------------------
[974]	582
	583
[1992]	584	! 2.1 - Undisturbed area and Wake integrals
	585	! ---------------------------------------------------------
[974]	586
[1992]	587	DO i = 1, klon
	588	z(i) = 0.
	589	ktop(i) = 0
	590	kupper(i) = 0
	591	sum_thu(i) = 0.
	592	sum_tu(i) = 0.
	593	sum_qu(i) = 0.
	594	sum_thvu(i) = 0.
	595	sum_dth(i) = 0.
	596	sum_dq(i) = 0.
	597	sum_rho(i) = 0.
	598	sum_dtdwn(i) = 0.
	599	sum_dqdwn(i) = 0.
[974]	600
[1992]	601	av_thu(i) = 0.
	602	av_tu(i) = 0.
	603	av_qu(i) = 0.
	604	av_thvu(i) = 0.
	605	av_dth(i) = 0.
	606	av_dq(i) = 0.
	607	av_rho(i) = 0.
	608	av_dtdwn(i) = 0.
	609	av_dqdwn(i) = 0.
	610	END DO
[974]	611
[1992]	612	! Distance between wakes
	613	DO i = 1, klon
	614	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
	615	END DO
	616	! Potential temperatures and humidity
	617	! ----------------------------------------------------------
	618	DO k = 1, klev
	619	DO i = 1, klon
	620	! write(,)'wake 1',i,k,rd,te(i,k)
	621	rho(i, k) = p(i, k)/(rd*te(i,k))
	622	! write(,)'wake 2',rho(i,k)
	623	IF (k==1) THEN
	624	! write(,)'wake 3',i,k,rd,te(i,k)
	625	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
	626	! write(,)'wake 4',i,k,rd,te(i,k)
	627	zhh(i, k) = 0
	628	ELSE
	629	! write(,)'wake 5',rd,(te(i,k)+te(i,k-1))
	630	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
	631	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
	632	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
	633	END IF
	634	! write(,)'wake 7',ppi(i,k)
	635	the(i, k) = te(i, k)/ppi(i, k)
	636	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
	637	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
	638	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
	639	! write(,)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
	640	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
	641	dth(i, k) = deltatw(i, k)/ppi(i, k)
	642	END DO
	643	END DO
[1403]	644
[1992]	645	DO k = 1, klev - 1
	646	DO i = 1, klon
	647	IF (k==1) THEN
	648	n2(i, k) = 0
	649	ELSE
[2635]	650	n2(i, k) = amax1(0., -rg*2/the(i,k)rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
	651	(p(i,k+1)-p(i,k-1)))
[1992]	652	END IF
	653	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
[1403]	654
[1992]	655	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
	656	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
	657	END DO
	658	END DO
[974]	659
[1992]	660	DO i = 1, klon
	661	n2(i, klev) = 0
	662	zh(i, klev) = 0
	663	cgw(i, klev) = 0
	664	tgw(i, klev) = 0
	665	END DO
[974]	666
[1992]	667	! Calcul de la masse volumique moyenne de la colonne (bdlmd)
	668	! -----------------------------------------------------------------
[974]	669
[1992]	670	DO k = 1, klev
	671	DO i = 1, klon
	672	epaisseur1(i, k) = 0.
	673	epaisseur2(i, k) = 0.
	674	END DO
	675	END DO
[974]	676
[1992]	677	DO i = 1, klon
	678	epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
	679	epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
	680	rhow_moyen(i, 1) = rhow(i, 1)
	681	END DO
[974]	682
[1992]	683	DO k = 2, klev
	684	DO i = 1, klon
	685	epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
	686	epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
	687	rhow_moyen(i, k) = (rhow_moyen(i,k-1)epaisseur2(i,k-1)+rhow(i,k) &
	688	epaisseur1(i,k))/epaisseur2(i, k)
	689	END DO
	690	END DO
[974]	691
	692
[1992]	693	! Choose an integration bound well above wake top
	694	! -----------------------------------------------------------------
[974]	695
[1992]	696	! Pupper = 50000. ! melting level
	697	! Pupper = 60000.
	698	! Pupper = 80000. ! essais pour case_e
	699	DO i = 1, klon
	700	pupper(i) = 0.6*ph(i, 1)
	701	pupper(i) = max(pupper(i), 45000.)
	702	! cc Pupper(i) = 60000.
	703	END DO
[1146]	704
[1403]	705
[1992]	706	! Determine Wake top pressure (Ptop) from buoyancy integral
	707	! --------------------------------------------------------
[1403]	708
[1992]	709	! -1/ Pressure of the level where dth becomes less than delta_t_min.
	710
	711	DO i = 1, klon
	712	ptop_provis(i) = ph(i, 1)
	713	END DO
	714	DO k = 2, klev
	715	DO i = 1, klon
	716
	717	! IM v3JYG; ptop_provis(i).LT. ph(i,1)
	718
	719	IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
	720	ptop_provis(i)==ph(i,1)) THEN
[2635]	721	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &
	722	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
[1992]	723	END IF
	724	END DO
	725	END DO
	726
	727	! -2/ dth integral
	728
	729	DO i = 1, klon
	730	sum_dth(i) = 0.
	731	dthmin(i) = -delta_t_min
	732	z(i) = 0.
	733	END DO
	734
	735	DO k = 1, klev
	736	DO i = 1, klon
	737	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
	738	IF (dz(i)>0) THEN
	739	z(i) = z(i) + dz(i)
	740	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	741	dthmin(i) = amin1(dthmin(i), dth(i,k))
	742	END IF
	743	END DO
	744	END DO
	745
	746	! -3/ height of triangle with area= sum_dth and base = dthmin
	747
	748	DO i = 1, klon
	749	hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
	750	hw0(i) = amax1(hwmin, hw0(i))
	751	END DO
	752
	753	! -4/ now, get Ptop
	754
	755	DO i = 1, klon
	756	z(i) = 0.
	757	ptop(i) = ph(i, 1)
	758	END DO
	759
	760	DO k = 1, klev
	761	DO i = 1, klon
	762	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
	763	IF (dz(i)>0) THEN
	764	z(i) = z(i) + dz(i)
	765	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
	766	END IF
	767	END DO
	768	END DO
	769
[2671]	770	IF (prt_level>=10) THEN
	771	PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
	772	ENDIF
[1992]	773
[2671]	774
[1992]	775	! -5/ Determination de ktop et kupper
	776
	777	DO k = klev, 1, -1
	778	DO i = 1, klon
	779	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	780	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
	781	END DO
	782	END DO
	783
	784	! On evite kupper = 1 et kupper = klev
	785	DO i = 1, klon
	786	kupper(i) = max(kupper(i), 2)
	787	kupper(i) = min(kupper(i), klev-1)
	788	END DO
	789
	790
	791	! -6/ Correct ktop and ptop
	792
	793	DO i = 1, klon
	794	ptop_new(i) = ptop(i)
	795	END DO
	796	DO k = klev, 2, -1
	797	DO i = 1, klon
	798	IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
	799	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
	800	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	801	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	802	END IF
	803	END DO
	804	END DO
	805
	806	DO i = 1, klon
	807	ptop(i) = ptop_new(i)
	808	END DO
	809
	810	DO k = klev, 1, -1
	811	DO i = 1, klon
	812	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	813	END DO
	814	END DO
	815
[2671]	816	IF (prt_level>=10) THEN
	817	PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
	818	ENDIF
	819
[1992]	820	! -5/ Set deltatw & deltaqw to 0 above kupper
	821
	822	DO k = 1, klev
	823	DO i = 1, klon
	824	IF (k>=kupper(i)) THEN
	825	deltatw(i, k) = 0.
	826	deltaqw(i, k) = 0.
[2635]	827	d_deltatw2(i,k) = -deltatw0(i,k)
	828	d_deltaqw2(i,k) = -deltaqw0(i,k)
[1992]	829	END IF
	830	END DO
	831	END DO
	832
	833
	834	! Vertical gradient of LS omega
	835
	836	DO k = 1, klev
	837	DO i = 1, klon
	838	IF (k<=kupper(i)) THEN
	839	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
	840	END IF
	841	END DO
	842	END DO
	843
	844	! Integrals (and wake top level number)
	845	! --------------------------------------
	846
	847	! Initialize sum_thvu to 1st level virt. pot. temp.
	848
	849	DO i = 1, klon
	850	z(i) = 1.
	851	dz(i) = 1.
[2495]	852	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
[1992]	853	sum_dth(i) = 0.
	854	END DO
	855
	856	DO k = 1, klev
	857	DO i = 1, klon
	858	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	859	IF (dz(i)>0) THEN
	860	z(i) = z(i) + dz(i)
	861	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	862	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	863	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
[2495]	864	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
[1992]	865	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	866	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	867	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	868	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	869	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	870	END IF
	871	END DO
	872	END DO
	873
	874	DO i = 1, klon
	875	hw0(i) = z(i)
	876	END DO
	877
	878
	879	! 2.1 - WAPE and mean forcing computation
	880	! ---------------------------------------
	881
	882	! ---------------------------------------
	883
	884	! Means
	885
	886	DO i = 1, klon
	887	av_thu(i) = sum_thu(i)/hw0(i)
	888	av_tu(i) = sum_tu(i)/hw0(i)
	889	av_qu(i) = sum_qu(i)/hw0(i)
	890	av_thvu(i) = sum_thvu(i)/hw0(i)
	891	! av_thve = sum_thve/hw0
	892	av_dth(i) = sum_dth(i)/hw0(i)
	893	av_dq(i) = sum_dq(i)/hw0(i)
	894	av_rho(i) = sum_rho(i)/hw0(i)
	895	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	896	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	897
[2635]	898	wape(i) = -rghw0(i)(av_dth(i)+ &
	899	epsim1(av_thu(i)av_dq(i)+av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
	900
[1992]	901	END DO
	902
	903	! 2.2 Prognostic variable update
	904	! ------------------------------
	905
	906	! Filter out bad wakes
	907
	908	DO k = 1, klev
	909	DO i = 1, klon
	910	IF (wape(i)<0.) THEN
	911	deltatw(i, k) = 0.
	912	deltaqw(i, k) = 0.
	913	dth(i, k) = 0.
[2635]	914	d_deltatw2(i,k) = -deltatw0(i,k)
	915	d_deltaqw2(i,k) = -deltaqw0(i,k)
[1992]	916	END IF
	917	END DO
	918	END DO
	919
	920	DO i = 1, klon
	921	IF (wape(i)<0.) THEN
	922	wape(i) = 0.
	923	cstar(i) = 0.
	924	hw(i) = hwmin
[2635]	925	!jyg<
	926	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
	927	sigmaw_targ = max(sigmad, sigd_con(i))
	928	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
	929	sigmaw(i) = sigmaw_targ
	930	!>jyg
[1992]	931	fip(i) = 0.
	932	gwake(i) = .FALSE.
	933	ELSE
[3208]	934	hw(i) = hw0(i)
[1992]	935	cstar(i) = starksqrt(2.wape(i))
	936	gwake(i) = .TRUE.
	937	END IF
	938	END DO
	939
	940
	941	! Check qx and qw positivity
	942	! --------------------------
	943	DO i = 1, klon
[2635]	944	q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), &
	945	(qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
[1992]	946	END DO
	947	DO k = 2, klev
	948	DO i = 1, klon
[2635]	949	q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
	950	(qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
[1992]	951	IF (q1_min(i)<=q0_min(i)) THEN
	952	q0_min(i) = q1_min(i)
	953	END IF
	954	END DO
	955	END DO
	956
	957	DO i = 1, klon
	958	ok_qx_qw(i) = q0_min(i) >= 0.
	959	alpha(i) = 1.
	960	END DO
	961
[2671]	962	IF (prt_level>=10) THEN
[2757]	963	PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
	964	sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
[2671]	965	ENDIF
	966
	967
[1992]	968	! C -----------------------------------------------------------------
	969	! Sub-time-stepping
	970	! -----------------
	971
	972	nsub = 10
	973	dtimesub = dtime/nsub
	974
	975	! ------------------------------------------------------------
	976	DO isubstep = 1, nsub
	977	! ------------------------------------------------------------
	978
	979	! wk_adv is the logical flag enabling wake evolution in the time advance
	980	! loop
	981	DO i = 1, klon
	982	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
	983	END DO
[2671]	984	IF (prt_level>=10) THEN
[2757]	985	PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
	986	isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
[2671]	987	ENDIF
[1992]	988
	989	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
	990	! négatif de ktop à kupper --------
	991	! cc On calcule pour cela une densité wdens0 pour laquelle on
	992	! aurait un entrainement nul ---
[3208]	993	!jyg<
	994	! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
	995	! les poches sont insuffisantes pour accueillir tout le flux de masse
	996	! des descentes unsaturees. Nous faisons alors l'hypothese que la
	997	! convection profonde cree directement de nouvelles poches, sans passer
	998	! par les thermiques. La nouvelle valeur de wdens est alors imposée.
	999
[1992]	1000	DO i = 1, klon
	1001	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
	1002	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
	1003	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
	1004	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
[2635]	1005	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
	1006	wdens0 = (sigmaw(i)/(4.3.14)) &
	1007	((1.-sigmaw(i))omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))cstar(i)))**(2)
[3252]	1008	IF (prt_level >= 10) THEN
	1009	print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
	1010	omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
	1011	ENDIF
[1992]	1012	IF (wdens(i)<=wdens0*1.1) THEN
[3208]	1013	IF (iflag_wk_pop_dyn >= 1) THEN
	1014	d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
	1015	ENDIF
[1992]	1016	wdens(i) = wdens0
	1017	END IF
	1018	END IF
	1019	END DO
	1020
	1021	DO i = 1, klon
	1022	IF (wk_adv(i)) THEN
[1403]	1023	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
[3208]	1024	rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
[2635]	1025	!jyg<
	1026	!! sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
	1027	sigmaw_targ = min(sigmaw(i), sigmaw_max)
	1028	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
	1029	sigmaw(i) = sigmaw_targ
	1030	!>jyg
[1992]	1031	END IF
	1032	END DO
[2635]	1033
[3208]	1034	IF (iflag_wk_pop_dyn >= 1) THEN
[3454]	1035	! The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
	1036	! Here, it has to be set to zero.
	1037	death_rate(:) = 0.
[3208]	1038
	1039	IF (iflag_wk_act ==2) THEN
	1040	DO i = 1, klon
	1041	IF (wk_adv(i)) THEN
	1042	wape1_act(i) = abs(cin(i))
	1043	wape2_act(i) = 2.*wape1_act(i) + 1.
	1044	act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
	1045	ENDIF ! (wk_adv(i))
	1046	ENDDO
	1047	ENDIF ! (iflag_wk_act ==2)
	1048
	1049
	1050	DO i = 1, klon
	1051	IF (wk_adv(i)) THEN
	1052	!! tau_wk(i) = max(rad_wk(i)/(3.cstar(i))((cstar(i)/cstart)**1.5 - 1), 100.)
	1053	tau_wk_inv(i) = max( (3.cstar(i))/(rad_wk(i)((cstar(i)/cstart)**1.5 - 1)), 0.)
	1054	tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
	1055	drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
[3252]	1056	(1 + 2f_shear(i)(2.sigmaw(i)-aa0wdens(i)) - 2.*sigmaw(i))
	1057	!! (1 - 2sigmaw(i)(1.-f_shear(i)))
[3208]	1058	drdt_pos=max(drdt(i),0.)
	1059
	1060	!! d_wdens(i) = ( wgen(i)(1.+2.(sigmaw(i)-sigmad)) &
	1061	!! - wdens(i)*tau_wk_inv_min &
	1062	!! - 2.gfl(i)wdens(i)Cstar(i) )dtimesub
	1063	d_awdens(i) = ( wgen(i) - (1./tau_cv)(awdens(i) - act(i)wdens(i)) )*dtimesub
	1064	d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min - &
	1065	2.wdens(i)gfl(i)drdt_pos )dtimesub
	1066	d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
	1067
	1068	!! d_sigmaw(i) = ( (1.-2f_shear(i)sigmaw(i))(gfl(i)Cstar(i)+wgen(i)*sigmad/wdens(i)) &
	1069	!! + 2.f_shear(i)wgen(i)sigmaw(i)*2/wdens(i) &
	1070	!! - sigmaw(i)tau_wk_inv_min )dtimesub
	1071	d_sig_gen(i) = wgen(i)*aa0
	1072	d_sig_death(i) = - sigmaw(i)(1.-awdens(i)/wdens(i))tau_wk_inv_min
	1073	!! d_sig_col(i) = - 2f_shear(i)sigmaw(i)gfl(i)drdt_pos
	1074	d_sig_col(i) = - 2f_shear(i)(2.sigmaw(i)-wdens(i)aa0)gfl(i)drdt_pos
	1075	d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)cstar(i) )dtimesub
	1076	d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
	1077	ENDIF
	1078	ENDDO
	1079
	1080	IF (prt_level >= 10) THEN
	1081	print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
	1082	cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
	1083	print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
	1084	wdens(1), awdens(1), act(1), d_awdens(1)
	1085	print ,'wake, wgen, -(wdens-awdens)tau_wk_inv, -2.wdensgfl*drdt_pos, d_wdens ', &
	1086	wgen(1), -(wdens(1)-awdens(1))tau_wk_inv(1), -2.wdens(1)gfl(1)drdt_pos, d_wdens(1)
	1087	print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
	1088	d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
	1089	ENDIF
	1090
	1091	ELSE ! (iflag_wk_pop_dyn >= 1)
	1092
	1093	! cc nrlmd
	1094
	1095	DO i = 1, klon
	1096	IF (wk_adv(i)) THEN
	1097	! cc nrlmd Introduction du taux de mortalité des poches et
	1098	! test sur sigmaw_max=0.4
	1099	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
	1100	IF (sigmaw(i)>=sigmaw_max) THEN
	1101	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
	1102	ELSE
	1103	death_rate(i) = 0.
	1104	END IF
	1105
	1106	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
	1107	dtimesub
	1108	! $ - nat_rate(i)sigmaw(i)dtimesub
	1109	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
	1110	! c $ death_rate(i),ktop(i),kupper(i)',
	1111	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
	1112	! c $ death_rate(i),ktop(i),kupper(i)
	1113
	1114	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
	1115	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
	1116	! wdens = wdens0/(10.*sigmaw)
	1117	! sigmaw =max(sigmaw,sigd_con)
	1118	! sigmaw =max(sigmaw,sigmad)
[1992]	1119	END IF
[3208]	1120	END DO
[2635]	1121
[3208]	1122	ENDIF ! (iflag_wk_pop_dyn >= 1)
[1403]	1123
[1992]	1124
	1125	! calcul de la difference de vitesse verticale poche - zone non perturbee
	1126	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
[2671]	1127	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
[1992]	1128	! IM 060208 au niveau k=1..?
[2671]	1129	!JYG 161013 Correction : maintenant omg est dimensionne a klev.
[1992]	1130	DO k = 1, klev
	1131	DO i = 1, klon
	1132	IF (wk_adv(i)) THEN !!! nrlmd
	1133	dp_deltomg(i, k) = 0.
	1134	END IF
	1135	END DO
	1136	END DO
[2671]	1137	DO k = 1, klev
[1992]	1138	DO i = 1, klon
	1139	IF (wk_adv(i)) THEN !!! nrlmd
	1140	omg(i, k) = 0.
	1141	END IF
	1142	END DO
	1143	END DO
	1144
	1145	DO i = 1, klon
	1146	IF (wk_adv(i)) THEN
	1147	z(i) = 0.
	1148	omg(i, 1) = 0.
	1149	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
	1150	END IF
	1151	END DO
	1152
	1153	DO k = 2, klev
	1154	DO i = 1, klon
	1155	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
[974]	1156	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
[1992]	1157	z(i) = z(i) + dz(i)
	1158	dp_deltomg(i, k) = dp_deltomg(i, 1)
	1159	omg(i, k) = dp_deltomg(i, 1)*z(i)
	1160	END IF
	1161	END DO
	1162	END DO
	1163
	1164	DO i = 1, klon
	1165	IF (wk_adv(i)) THEN
	1166	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
	1167	ztop(i) = z(i) + dztop(i)
	1168	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
	1169	END IF
	1170	END DO
	1171
[2671]	1172	IF (prt_level>=10) THEN
	1173	PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
[2757]	1174	PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
	1175	omgtop(igout), ptop(igout), ktop(igout)
[2671]	1176	ENDIF
	1177
[1992]	1178	! -----------------
	1179	! From m/s to Pa/s
	1180	! -----------------
	1181
	1182	DO i = 1, klon
	1183	IF (wk_adv(i)) THEN
	1184	omgtop(i) = -rho(i, ktop(i))rgomgtop(i)
	1185	dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
	1186	END IF
	1187	END DO
	1188
	1189	DO k = 1, klev
	1190	DO i = 1, klon
	1191	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
	1192	omg(i, k) = -rho(i, k)rgomg(i, k)
	1193	dp_deltomg(i, k) = dp_deltomg(i, 1)
	1194	END IF
	1195	END DO
	1196	END DO
	1197
	1198	! raccordement lineaire de omg de ptop a pupper
	1199
	1200	DO i = 1, klon
	1201	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
	1202	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
	1203	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
	1204	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
	1205	(ptop(i)-pupper(i))
	1206	END IF
	1207	END DO
	1208
	1209	! c DO i=1,klon
	1210	! c print*,'Pente entre 0 et kupper (référence)'
	1211	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
	1212	! c print*,'Pente entre ktop et kupper'
	1213	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
	1214	! c ENDDO
	1215	! c
	1216	DO k = 1, klev
	1217	DO i = 1, klon
	1218	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
	1219	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
	1220	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
	1221	END IF
	1222	END DO
	1223	END DO
[2671]	1224	!! print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
[1992]	1225	! cc nrlmd
	1226	! c DO i=1,klon
	1227	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
	1228	! c END DO
	1229	! cc
	1230
	1231
	1232	! -- Compute wake average vertical velocity omgbw
	1233
	1234
[2671]	1235	DO k = 1, klev
[1992]	1236	DO i = 1, klon
[1146]	1237	IF (wk_adv(i)) THEN
[1992]	1238	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
	1239	END IF
	1240	END DO
	1241	END DO
	1242	! -- and its vertical gradient dp_omgbw
	1243
[2671]	1244	DO k = 1, klev-1
[1992]	1245	DO i = 1, klon
[1146]	1246	IF (wk_adv(i)) THEN
[1992]	1247	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
	1248	END IF
	1249	END DO
	1250	END DO
[2671]	1251	DO i = 1, klon
	1252	IF (wk_adv(i)) THEN
	1253	dp_omgbw(i, klev) = 0.
	1254	END IF
	1255	END DO
[974]	1256
[1992]	1257	! -- Upstream coefficients for omgb velocity
	1258	! -- (alpha_up(k) is the coefficient of the value at level k)
	1259	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
	1260	DO k = 1, klev
	1261	DO i = 1, klon
	1262	IF (wk_adv(i)) THEN
	1263	alpha_up(i, k) = 0.
	1264	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
	1265	END IF
	1266	END DO
	1267	END DO
[974]	1268
[1992]	1269	! Matrix expressing [The,deltatw] from [Th1,Th2]
[974]	1270
[1992]	1271	DO i = 1, klon
	1272	IF (wk_adv(i)) THEN
	1273	rre1(i) = 1. - sigmaw(i)
	1274	rre2(i) = sigmaw(i)
	1275	END IF
	1276	END DO
	1277	rrd1 = -1.
	1278	rrd2 = 1.
[974]	1279
[1992]	1280	! -- Get [Th1,Th2], dth and [q1,q2]
[974]	1281
[1992]	1282	DO k = 1, klev
	1283	DO i = 1, klon
	1284	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
	1285	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1286	th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
	1287	th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
	1288	q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
	1289	q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
	1290	END IF
	1291	END DO
	1292	END DO
[974]	1293
[1992]	1294	DO i = 1, klon
	1295	IF (wk_adv(i)) THEN !!! nrlmd
	1296	d_th1(i, 1) = 0.
	1297	d_th2(i, 1) = 0.
	1298	d_dth(i, 1) = 0.
	1299	d_q1(i, 1) = 0.
	1300	d_q2(i, 1) = 0.
	1301	d_dq(i, 1) = 0.
	1302	END IF
	1303	END DO
[974]	1304
[1992]	1305	DO k = 2, klev
	1306	DO i = 1, klon
	1307	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
	1308	d_th1(i, k) = th1(i, k-1) - th1(i, k)
	1309	d_th2(i, k) = th2(i, k-1) - th2(i, k)
	1310	d_dth(i, k) = dth(i, k-1) - dth(i, k)
	1311	d_q1(i, k) = q1(i, k-1) - q1(i, k)
	1312	d_q2(i, k) = q2(i, k-1) - q2(i, k)
	1313	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
	1314	END IF
	1315	END DO
	1316	END DO
[1146]	1317
[1992]	1318	DO i = 1, klon
	1319	IF (wk_adv(i)) THEN
	1320	omgbdth(i, 1) = 0.
	1321	omgbdq(i, 1) = 0.
	1322	END IF
	1323	END DO
[1277]	1324
[1992]	1325	DO k = 2, klev
	1326	DO i = 1, klon
	1327	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
	1328	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
	1329	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
	1330	END IF
	1331	END DO
	1332	END DO
[1403]	1333
[2671]	1334	IF (prt_level>=10) THEN
	1335	PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev)
	1336	PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev)
	1337	PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev)
	1338	PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev)
	1339	ENDIF
	1340
[1992]	1341	! -----------------------------------------------------------------
[2671]	1342	DO k = 1, klev-1
[1992]	1343	DO i = 1, klon
	1344	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
	1345	! -----------------------------------------------------------------
[974]	1346
[1992]	1347	! Compute redistribution (advective) term
[1403]	1348
[1992]	1349	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
[2635]	1350	(rrd1omg(i,k)sigmaw(i)*d_th1(i,k) - &
	1351	rrd2omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1)- &
	1352	(1.-alpha_up(i,k))*omgbdth(i,k)- &
	1353	alpha_up(i,k+1)omgbdth(i,k+1))ppi(i, k)
[2671]	1354	! print*,'d_deltatw=', k, d_deltatw(i,k)
[1403]	1355
[1992]	1356	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
[2635]	1357	(rrd1omg(i,k)sigmaw(i)*d_q1(i,k)- &
	1358	rrd2omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1)- &
	1359	(1.-alpha_up(i,k))*omgbdq(i,k)- &
	1360	alpha_up(i,k+1)*omgbdq(i,k+1))
[2671]	1361	! print*,'d_deltaqw=', k, d_deltaqw(i,k)
[974]	1362
[1992]	1363	! and increment large scale tendencies
[974]	1364
	1365
	1366
	1367
[1992]	1368	! C
	1369	! -----------------------------------------------------------------
[2635]	1370	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i,k)- &
	1371	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/ &
	1372	(ph(i,k)-ph(i,k+1)) &
	1373	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
	1374	(ph(i,k)-ph(i,k+1)) )*ppi(i, k)
[974]	1375
[2635]	1376	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)- &
	1377	rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/ &
	1378	(ph(i,k)-ph(i,k+1)) &
	1379	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
	1380	(ph(i,k)-ph(i,k+1)) )
[1992]	1381	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
[2635]	1382	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]	1383
[2635]	1384	d_qe(i, k) = dtimesub(rre1(i)omg(i,k)sigmaw(i)d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
[1403]	1385
[1992]	1386	END IF
	1387	! cc
	1388	END DO
	1389	END DO
	1390	! ------------------------------------------------------------------
[974]	1391
[2671]	1392	IF (prt_level>=10) THEN
	1393	PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
	1394	PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
	1395	ENDIF
	1396
[1992]	1397	! Increment state variables
[3208]	1398	!jyg<
	1399	IF (iflag_wk_pop_dyn >= 1) THEN
	1400	DO k = 1, klev
	1401	DO i = 1, klon
	1402	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1403	detr(i,k) = - d_sig_death(i) - d_sig_col(i)
	1404	entr(i,k) = d_sig_gen(i)
	1405	ENDIF
	1406	ENDDO
	1407	ENDDO
	1408	ELSE ! (iflag_wk_pop_dyn >= 1)
	1409	DO k = 1, klev
	1410	DO i = 1, klon
	1411	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1412	detr(i, k) = 0.
	1413
	1414	entr(i, k) = 0.
	1415	ENDIF
	1416	ENDDO
	1417	ENDDO
	1418	ENDIF ! (iflag_wk_pop_dyn >= 1)
[974]	1419
[3208]	1420
	1421
[1992]	1422	DO k = 1, klev
	1423	DO i = 1, klon
	1424	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
	1425	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1426	! cc
[974]	1427
[1146]	1428
[974]	1429
[1992]	1430	! Coefficient de répartition
[974]	1431
[1992]	1432	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
	1433	(ph(i,kupper(i))-ph(i,1))
[2635]	1434	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
	1435	(p(i,1)-ph(i,kupper(i)))
[974]	1436
	1437
[1992]	1438	! Reintroduce compensating subsidence term.
[1146]	1439
[1992]	1440	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
	1441	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
	1442	! . /(1-sigmaw)
	1443	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
	1444	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
	1445	! . /(1-sigmaw)
[974]	1446
[1992]	1447	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
	1448	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
	1449	! . /(1-sigmaw)
	1450	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
	1451	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
	1452	! . /(1-sigmaw)
[974]	1453
[1992]	1454	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
	1455	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
	1456	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
[974]	1457
[2155]	1458	!
[1146]	1459
[1992]	1460	! cc nrlmd Prise en compte du taux de mortalité
	1461	! cc Définitions de entr, detr
[3208]	1462	!jyg<
	1463	!! detr(i, k) = 0.
	1464	!!
	1465	!! entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
	1466	!! sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
	1467	!!
	1468	entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
	1469	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
	1470	!>jyg
	1471	spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
[1146]	1472
[1992]	1473	! cc spread(i,k) =
	1474	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
	1475	! cc $ sigmaw(i)
[1146]	1476
	1477
[1992]	1478	! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
	1479	! Jingmei
[1146]	1480
[1992]	1481	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
	1482	! & Tgw(i,k),deltatw(i,k)
	1483	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
	1484	dtimesub
	1485	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
	1486	ff(i) = d_deltatw(i, k)/dtimesub
[1403]	1487
[1992]	1488	! Sans GW
[1403]	1489
[1992]	1490	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
[974]	1491
[1992]	1492	! GW formule 1
	1493
	1494	! deltatw(k) = deltatw(k)+dtimesub*
	1495	! $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
	1496
	1497	! GW formule 2
	1498
	1499	IF (dtimesub*tgw(i,k)<1.E-10) THEN
[2635]	1500	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
	1501	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
	1502	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
	1503	tgw(i,k)*deltatw(i,k) )
[1992]	1504	ELSE
[2635]	1505	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i,k)))* &
	1506	(ff(i)+dtke(i,k) - &
	1507	entr(i,k)*deltatw(i,k)/sigmaw(i) - &
	1508	(death_rate(i)sigmaw(i)+detr(i,k))deltatw(i,k)/(1.-sigmaw(i)) - &
	1509	tgw(i,k)*deltatw(i,k) )
[1992]	1510	END IF
	1511
	1512	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1513
	1514	gg(i) = d_deltaqw(i, k)/dtimesub
	1515
[2635]	1516	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
	1517	entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
	1518	(death_rate(i)sigmaw(i)+detr(i,k))deltaqw(i,k)/(1.-sigmaw(i)))
[1992]	1519	! cc
	1520
	1521	! cc nrlmd
	1522	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
	1523	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
	1524	! cc
	1525	END IF
	1526	END DO
	1527	END DO
	1528
	1529
	1530	! Scale tendencies so that water vapour remains positive in w and x.
	1531
	1532	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
	1533	d_deltaqw, sigmaw, d_sigmaw, alpha)
	1534
	1535	! cc nrlmd
	1536	! c print*,'alpha'
	1537	! c do i=1,klon
	1538	! c print*,alpha(i)
	1539	! c end do
	1540	! cc
	1541	DO k = 1, klev
	1542	DO i = 1, klon
	1543	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1544	d_te(i, k) = alpha(i)*d_te(i, k)
	1545	d_qe(i, k) = alpha(i)*d_qe(i, k)
	1546	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
	1547	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
	1548	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
	1549	END IF
	1550	END DO
	1551	END DO
	1552	DO i = 1, klon
	1553	IF (wk_adv(i)) THEN
	1554	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
	1555	END IF
	1556	END DO
	1557
	1558	! Update large scale variables and wake variables
	1559	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
	1560	! IM 060208 DO k = 1,kupper(i)
	1561	DO k = 1, klev
	1562	DO i = 1, klon
	1563	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1564	dtls(i, k) = dtls(i, k) + d_te(i, k)
	1565	dqls(i, k) = dqls(i, k) + d_qe(i, k)
	1566	! cc nrlmd
	1567	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
	1568	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
	1569	! cc
	1570	END IF
	1571	END DO
	1572	END DO
	1573	DO k = 1, klev
	1574	DO i = 1, klon
	1575	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1576	te(i, k) = te0(i, k) + dtls(i, k)
	1577	qe(i, k) = qe0(i, k) + dqls(i, k)
	1578	the(i, k) = te(i, k)/ppi(i, k)
	1579	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
	1580	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
	1581	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1582	! c print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
	1583	! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
	1584	END IF
	1585	END DO
	1586	END DO
[3208]	1587	!
[1992]	1588	DO i = 1, klon
	1589	IF (wk_adv(i)) THEN
	1590	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
[2635]	1591	d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
[1992]	1592	END IF
	1593	END DO
[3208]	1594	!jyg<
	1595	IF (iflag_wk_pop_dyn >= 1) THEN
	1596	DO i = 1, klon
	1597	IF (wk_adv(i)) THEN
	1598	awdens(i) = awdens(i) + d_awdens(i)
	1599	wdens(i) = wdens(i) + d_wdens(i)
	1600	d_awdens2(i) = d_awdens2(i) + d_awdens(i)
	1601	d_wdens2(i) = d_wdens2(i) + d_wdens(i)
	1602	END IF
	1603	END DO
	1604	DO i = 1, klon
	1605	IF (wk_adv(i)) THEN
	1606	wdens_targ = max(wdens(i),wdensmin)
	1607	d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
	1608	wdens(i) = wdens_targ
	1609	!
	1610	wdens_targ = min( max(awdens(i),0.), wdens(i) )
	1611	d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
	1612	awdens(i) = wdens_targ
	1613	END IF
	1614	END DO
	1615	DO i = 1, klon
	1616	IF (wk_adv(i)) THEN
	1617	sigmaw_targ = max(sigmaw(i),sigmad)
	1618	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
	1619	sigmaw(i) = sigmaw_targ
	1620	END IF
	1621	END DO
	1622	ENDIF ! (iflag_wk_pop_dyn >= 1)
	1623	!>jyg
[1992]	1624
	1625
	1626	! Determine Ptop from buoyancy integral
	1627	! ---------------------------------------
	1628
	1629	! - 1/ Pressure of the level where dth changes sign.
	1630
	1631	DO i = 1, klon
	1632	IF (wk_adv(i)) THEN
	1633	ptop_provis(i) = ph(i, 1)
	1634	END IF
	1635	END DO
	1636
	1637	DO k = 2, klev
	1638	DO i = 1, klon
	1639	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
	1640	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
[2635]	1641	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
	1642	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
[1992]	1643	END IF
	1644	END DO
	1645	END DO
	1646
	1647	! - 2/ dth integral
	1648
	1649	DO i = 1, klon
	1650	IF (wk_adv(i)) THEN !!! nrlmd
	1651	sum_dth(i) = 0.
	1652	dthmin(i) = -delta_t_min
[974]	1653	z(i) = 0.
[1992]	1654	END IF
	1655	END DO
	1656
	1657	DO k = 1, klev
	1658	DO i = 1, klon
	1659	IF (wk_adv(i)) THEN
	1660	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
	1661	IF (dz(i)>0) THEN
	1662	z(i) = z(i) + dz(i)
	1663	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1664	dthmin(i) = amin1(dthmin(i), dth(i,k))
	1665	END IF
	1666	END IF
	1667	END DO
	1668	END DO
	1669
	1670	! - 3/ height of triangle with area= sum_dth and base = dthmin
	1671
	1672	DO i = 1, klon
	1673	IF (wk_adv(i)) THEN
	1674	hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
	1675	hw(i) = amax1(hwmin, hw(i))
	1676	END IF
	1677	END DO
	1678
	1679	! - 4/ now, get Ptop
	1680
	1681	DO i = 1, klon
	1682	IF (wk_adv(i)) THEN !!! nrlmd
	1683	ktop(i) = 0
	1684	z(i) = 0.
	1685	END IF
	1686	END DO
	1687
	1688	DO k = 1, klev
	1689	DO i = 1, klon
	1690	IF (wk_adv(i)) THEN
	1691	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
	1692	IF (dz(i)>0) THEN
	1693	z(i) = z(i) + dz(i)
	1694	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
	1695	ktop(i) = k
	1696	END IF
	1697	END IF
	1698	END DO
	1699	END DO
	1700
	1701	! 4.5/Correct ktop and ptop
	1702
	1703	DO i = 1, klon
	1704	IF (wk_adv(i)) THEN
	1705	ptop_new(i) = ptop(i)
	1706	END IF
	1707	END DO
	1708
	1709	DO k = klev, 2, -1
	1710	DO i = 1, klon
	1711	! IM v3JYG; IF (k .GE. ktop(i)
	1712	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
	1713	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
[2635]	1714	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
	1715	(dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
[1992]	1716	END IF
	1717	END DO
	1718	END DO
	1719
	1720
	1721	DO i = 1, klon
	1722	IF (wk_adv(i)) THEN
	1723	ptop(i) = ptop_new(i)
	1724	END IF
	1725	END DO
	1726
	1727	DO k = klev, 1, -1
	1728	DO i = 1, klon
	1729	IF (wk_adv(i)) THEN !!! nrlmd
	1730	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	1731	END IF
	1732	END DO
	1733	END DO
	1734
	1735	! 5/ Set deltatw & deltaqw to 0 above kupper
	1736
	1737	DO k = 1, klev
	1738	DO i = 1, klon
	1739	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
	1740	deltatw(i, k) = 0.
	1741	deltaqw(i, k) = 0.
[2635]	1742	d_deltatw2(i,k) = -deltatw0(i,k)
	1743	d_deltaqw2(i,k) = -deltaqw0(i,k)
[1992]	1744	END IF
	1745	END DO
	1746	END DO
	1747
	1748
	1749	! -------------Cstar computation---------------------------------
	1750	DO i = 1, klon
	1751	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1752	sum_thu(i) = 0.
	1753	sum_tu(i) = 0.
	1754	sum_qu(i) = 0.
	1755	sum_thvu(i) = 0.
	1756	sum_dth(i) = 0.
	1757	sum_dq(i) = 0.
	1758	sum_rho(i) = 0.
	1759	sum_dtdwn(i) = 0.
	1760	sum_dqdwn(i) = 0.
	1761
	1762	av_thu(i) = 0.
[1992]	1763	av_tu(i) = 0.
	1764	av_qu(i) = 0.
[974]	1765	av_thvu(i) = 0.
	1766	av_dth(i) = 0.
	1767	av_dq(i) = 0.
[1992]	1768	av_rho(i) = 0.
	1769	av_dtdwn(i) = 0.
[974]	1770	av_dqdwn(i) = 0.
[1992]	1771	END IF
	1772	END DO
[974]	1773
[1992]	1774	! Integrals (and wake top level number)
	1775	! --------------------------------------
[974]	1776
[1992]	1777	! Initialize sum_thvu to 1st level virt. pot. temp.
[974]	1778
[1992]	1779	DO i = 1, klon
	1780	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1781	z(i) = 1.
	1782	dz(i) = 1.
[2495]	1783	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
[974]	1784	sum_dth(i) = 0.
[1992]	1785	END IF
	1786	END DO
[974]	1787
[1992]	1788	DO k = 1, klev
	1789	DO i = 1, klon
	1790	IF (wk_adv(i)) THEN !!! nrlmd
	1791	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	1792	IF (dz(i)>0) THEN
	1793	z(i) = z(i) + dz(i)
	1794	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	1795	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	1796	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
[2495]	1797	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
[1992]	1798	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1799	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	1800	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	1801	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	1802	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	1803	END IF
	1804	END IF
	1805	END DO
	1806	END DO
	1807
	1808	DO i = 1, klon
	1809	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1810	hw0(i) = z(i)
[1992]	1811	END IF
	1812	END DO
[974]	1813
	1814
[1992]	1815	! - WAPE and mean forcing computation
	1816	! ---------------------------------------
	1817
	1818	! ---------------------------------------
	1819
	1820	! Means
	1821
	1822	DO i = 1, klon
	1823	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1824	av_thu(i) = sum_thu(i)/hw0(i)
	1825	av_tu(i) = sum_tu(i)/hw0(i)
	1826	av_qu(i) = sum_qu(i)/hw0(i)
	1827	av_thvu(i) = sum_thvu(i)/hw0(i)
	1828	av_dth(i) = sum_dth(i)/hw0(i)
	1829	av_dq(i) = sum_dq(i)/hw0(i)
	1830	av_rho(i) = sum_rho(i)/hw0(i)
	1831	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	1832	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	1833
[2635]	1834	wape(i) = -rghw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
	1835	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
[1992]	1836	END IF
	1837	END DO
[974]	1838
[1992]	1839	! Filter out bad wakes
[974]	1840
[1992]	1841	DO k = 1, klev
	1842	DO i = 1, klon
	1843	IF (wk_adv(i)) THEN !!! nrlmd
	1844	IF (wape(i)<0.) THEN
	1845	deltatw(i, k) = 0.
	1846	deltaqw(i, k) = 0.
	1847	dth(i, k) = 0.
[2635]	1848	d_deltatw2(i,k) = -deltatw0(i,k)
	1849	d_deltaqw2(i,k) = -deltaqw0(i,k)
[1992]	1850	END IF
	1851	END IF
	1852	END DO
	1853	END DO
[974]	1854
[1992]	1855	DO i = 1, klon
	1856	IF (wk_adv(i)) THEN !!! nrlmd
	1857	IF (wape(i)<0.) THEN
	1858	wape(i) = 0.
	1859	cstar(i) = 0.
	1860	hw(i) = hwmin
[2635]	1861	!jyg<
	1862	!! sigmaw(i) = max(sigmad, sigd_con(i))
	1863	sigmaw_targ = max(sigmad, sigd_con(i))
	1864	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
	1865	sigmaw(i) = sigmaw_targ
	1866	!>jyg
[1992]	1867	fip(i) = 0.
	1868	gwake(i) = .FALSE.
	1869	ELSE
	1870	cstar(i) = starksqrt(2.wape(i))
	1871	gwake(i) = .TRUE.
	1872	END IF
	1873	END IF
	1874	END DO
	1875
	1876	END DO ! end sub-timestep loop
	1877
[2671]	1878	IF (prt_level>=10) THEN
[2757]	1879	PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
	1880	sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
[2671]	1881	ENDIF
[1992]	1882
	1883
	1884	! ----------------------------------------------------------
	1885	! Determine wake final state; recompute wape, cstar, ktop;
	1886	! filter out bad wakes.
	1887	! ----------------------------------------------------------
	1888
	1889	! 2.1 - Undisturbed area and Wake integrals
	1890	! ---------------------------------------------------------
	1891
	1892	DO i = 1, klon
	1893	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
	1894	IF (ok_qx_qw(i)) THEN
	1895	! cc
	1896	z(i) = 0.
	1897	sum_thu(i) = 0.
	1898	sum_tu(i) = 0.
	1899	sum_qu(i) = 0.
	1900	sum_thvu(i) = 0.
	1901	sum_dth(i) = 0.
[2757]	1902	sum_half_dth(i) = 0.
[1992]	1903	sum_dq(i) = 0.
	1904	sum_rho(i) = 0.
	1905	sum_dtdwn(i) = 0.
	1906	sum_dqdwn(i) = 0.
	1907
	1908	av_thu(i) = 0.
	1909	av_tu(i) = 0.
	1910	av_qu(i) = 0.
	1911	av_thvu(i) = 0.
	1912	av_dth(i) = 0.
	1913	av_dq(i) = 0.
	1914	av_rho(i) = 0.
	1915	av_dtdwn(i) = 0.
	1916	av_dqdwn(i) = 0.
[2757]	1917
	1918	dthmin(i) = -delta_t_min
[1992]	1919	END IF
	1920	END DO
	1921	! Potential temperatures and humidity
	1922	! ----------------------------------------------------------
	1923
	1924	DO k = 1, klev
	1925	DO i = 1, klon
	1926	! cc nrlmd IF ( wk_adv(i)) THEN
	1927	IF (ok_qx_qw(i)) THEN
	1928	! cc
	1929	rho(i, k) = p(i, k)/(rd*te(i,k))
	1930	IF (k==1) THEN
	1931	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
	1932	zhh(i, k) = 0
	1933	ELSE
	1934	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
	1935	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
	1936	END IF
	1937	the(i, k) = te(i, k)/ppi(i, k)
	1938	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
	1939	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
	1940	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
	1941	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
	1942	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1943	END IF
	1944	END DO
	1945	END DO
	1946
	1947	! Integrals (and wake top level number)
	1948	! -----------------------------------------------------------
	1949
	1950	! Initialize sum_thvu to 1st level virt. pot. temp.
	1951
	1952	DO i = 1, klon
	1953	! cc nrlmd IF ( wk_adv(i)) THEN
	1954	IF (ok_qx_qw(i)) THEN
	1955	! cc
	1956	z(i) = 1.
	1957	dz(i) = 1.
[2757]	1958	dz_half(i) = 1.
[2495]	1959	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
[1992]	1960	sum_dth(i) = 0.
	1961	END IF
	1962	END DO
	1963
	1964	DO k = 1, klev
	1965	DO i = 1, klon
	1966	! cc nrlmd IF ( wk_adv(i)) THEN
	1967	IF (ok_qx_qw(i)) THEN
	1968	! cc
	1969	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
[2757]	1970	dz_half(i) = -(amax1(ph(i,k+1),0.5(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)rg)
[1992]	1971	IF (dz(i)>0) THEN
	1972	z(i) = z(i) + dz(i)
	1973	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	1974	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	1975	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
[2495]	1976	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
[1992]	1977	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1978	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	1979	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	1980	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	1981	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
[2757]	1982	!
	1983	dthmin(i) = min(dthmin(i), dth(i,k))
[1992]	1984	END IF
[2757]	1985	IF (dz_half(i)>0) THEN
	1986	sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
	1987	END IF
[1992]	1988	END IF
	1989	END DO
	1990	END DO
	1991
	1992	DO i = 1, klon
	1993	! cc nrlmd IF ( wk_adv(i)) THEN
	1994	IF (ok_qx_qw(i)) THEN
	1995	! cc
	1996	hw0(i) = z(i)
	1997	END IF
	1998	END DO
	1999
	2000	! - WAPE and mean forcing computation
	2001	! -------------------------------------------------------------
	2002
	2003	! Means
	2004
	2005	DO i = 1, klon
	2006	! cc nrlmd IF ( wk_adv(i)) THEN
	2007	IF (ok_qx_qw(i)) THEN
	2008	! cc
	2009	av_thu(i) = sum_thu(i)/hw0(i)
	2010	av_tu(i) = sum_tu(i)/hw0(i)
	2011	av_qu(i) = sum_qu(i)/hw0(i)
	2012	av_thvu(i) = sum_thvu(i)/hw0(i)
	2013	av_dth(i) = sum_dth(i)/hw0(i)
	2014	av_dq(i) = sum_dq(i)/hw0(i)
	2015	av_rho(i) = sum_rho(i)/hw0(i)
	2016	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	2017	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	2018
[2635]	2019	wape2(i) = -rghw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i) + &
	2020	av_dth(i)av_qu(i)+av_dth(i)av_dq(i)))/av_thvu(i)
[1992]	2021	END IF
	2022	END DO
	2023
[2635]	2024
	2025
[1992]	2026	! Prognostic variable update
	2027	! ------------------------------------------------------------
	2028
	2029	! Filter out bad wakes
	2030
[2922]	2031	IF (iflag_wk_check_trgl>=1) THEN
[2757]	2032	! Check triangular shape of dth profile
	2033	DO i = 1, klon
	2034	IF (ok_qx_qw(i)) THEN
	2035	!! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
	2036	!! print ,'wake, 2.sum_dth(i)/(hw0(i)*dthmin(i)) ', &
	2037	!! 2.sum_dth(i)/(hw0(i)dthmin(i))
	2038	!! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
	2039	!! sum_half_dth(i), sum_dth(i)
	2040	IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
	2041	wape2(i) = -1.
	2042	!! print *,'wake, rej 1'
[2922]	2043	ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.sum_dth(i)/(hw0(i)dthmin(i)) - 1.) > 0.5) THEN
[2757]	2044	wape2(i) = -1.
	2045	!! print *,'wake, rej 2'
	2046	ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
	2047	wape2(i) = -1.
	2048	!! print *,'wake, rej 3'
	2049	END IF
	2050	END IF
	2051	END DO
	2052	END IF
	2053
	2054
[1992]	2055	DO k = 1, klev
	2056	DO i = 1, klon
	2057	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
	2058	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
	2059	! cc
	2060	deltatw(i, k) = 0.
	2061	deltaqw(i, k) = 0.
	2062	dth(i, k) = 0.
[2635]	2063	d_deltatw2(i,k) = -deltatw0(i,k)
	2064	d_deltaqw2(i,k) = -deltaqw0(i,k)
[1992]	2065	END IF
	2066	END DO
	2067	END DO
	2068
	2069
	2070	DO i = 1, klon
	2071	! cc nrlmd IF ( wk_adv(i)) THEN
	2072	IF (ok_qx_qw(i)) THEN
	2073	! cc
	2074	IF (wape2(i)<0.) THEN
[974]	2075	wape2(i) = 0.
[1992]	2076	cstar2(i) = 0.
[974]	2077	hw(i) = hwmin
[2635]	2078	!jyg<
	2079	!! sigmaw(i) = amax1(sigmad, sigd_con(i))
	2080	sigmaw_targ = max(sigmad, sigd_con(i))
	2081	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
	2082	sigmaw(i) = sigmaw_targ
	2083	!>jyg
[974]	2084	fip(i) = 0.
	2085	gwake(i) = .FALSE.
	2086	ELSE
[1992]	2087	IF (prt_level>=10) PRINT *, 'wape2>0'
	2088	cstar2(i) = starksqrt(2.wape2(i))
[974]	2089	gwake(i) = .TRUE.
[1992]	2090	END IF
	2091	END IF
	2092	END DO
[974]	2093
[1992]	2094	DO i = 1, klon
	2095	! cc nrlmd IF ( wk_adv(i)) THEN
	2096	IF (ok_qx_qw(i)) THEN
	2097	! cc
	2098	ktopw(i) = ktop(i)
	2099	END IF
	2100	END DO
[974]	2101
[1992]	2102	DO i = 1, klon
	2103	! cc nrlmd IF ( wk_adv(i)) THEN
	2104	IF (ok_qx_qw(i)) THEN
	2105	! cc
	2106	IF (ktopw(i)>0 .AND. gwake(i)) THEN
[1403]	2107
[1992]	2108	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
	2109	! cc heff = 600.
	2110	! Utilisation de la hauteur hw
	2111	! c heff = 0.7*hw
	2112	heff(i) = hw(i)
[1403]	2113
[1992]	2114	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
	2115	sqrt(sigmaw(i)wdens(i)3.14)
	2116	fip(i) = alpk*fip(i)
	2117	! jyg2
	2118	ELSE
	2119	fip(i) = 0.
	2120	END IF
	2121	END IF
	2122	END DO
[1146]	2123
[1992]	2124	! Limitation de sigmaw
	2125
	2126	! cc nrlmd
	2127	! DO i=1,klon
	2128	! IF (OK_qx_qw(i)) THEN
	2129	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
	2130	! ENDIF
	2131	! ENDDO
	2132	! cc
[3208]	2133
	2134	!jyg<
	2135	IF (iflag_wk_pop_dyn >= 1) THEN
	2136	DO i = 1, klon
	2137	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
	2138	.NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
	2139	ENDDO
	2140	ELSE ! (iflag_wk_pop_dyn >= 1)
	2141	DO i = 1, klon
	2142	kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
	2143	.NOT. ok_qx_qw(i)
	2144	ENDDO
	2145	ENDIF ! (iflag_wk_pop_dyn >= 1)
	2146	!>jyg
	2147
[1992]	2148	DO k = 1, klev
	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
[1992]	2153	! cc
	2154	dtls(i, k) = 0.
	2155	dqls(i, k) = 0.
	2156	deltatw(i, k) = 0.
	2157	deltaqw(i, k) = 0.
[2635]	2158	d_deltatw2(i,k) = -deltatw0(i,k)
	2159	d_deltaqw2(i,k) = -deltaqw0(i,k)
[3208]	2160	END IF ! (kill_wake(i))
[1992]	2161	END DO
	2162	END DO
	2163
	2164	DO i = 1, klon
[3208]	2165	!!jyg IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
	2166	!!jyg .NOT. ok_qx_qw(i)) THEN
	2167	IF (kill_wake(i)) THEN
[2635]	2168	ktopw(i) = 0
[1992]	2169	wape(i) = 0.
	2170	cstar(i) = 0.
[3208]	2171	!!jyg Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
[2308]	2172	!! hw(i) = hwmin !jyg
	2173	!! sigmaw(i) = sigmad !jyg
	2174	hw(i) = 0. !jyg
[1992]	2175	fip(i) = 0.
[3208]	2176	!! sigmaw(i) = 0. !jyg
	2177	sigmaw_targ = 0.
	2178	d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
	2179	sigmaw(i) = sigmaw_targ
	2180	IF (iflag_wk_pop_dyn >= 1) THEN
	2181	!! awdens(i) = 0.
	2182	!! wdens(i) = 0.
	2183	wdens_targ = 0.
	2184	d_wdens2(i) = wdens_targ - wdens(i)
	2185	wdens(i) = wdens_targ
	2186	wdens_targ = 0.
	2187	d_awdens2(i) = wdens_targ - awdens(i)
	2188	awdens(i) = wdens_targ
	2189	ENDIF ! (iflag_wk_pop_dyn >= 1)
	2190	ELSE ! (kill_wake(i))
[1992]	2191	wape(i) = wape2(i)
	2192	cstar(i) = cstar2(i)
[3208]	2193	END IF ! (kill_wake(i))
[1992]	2194	! c print*,'wape wape2 ktopw OK_qx_qw =',
	2195	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
	2196	END DO
	2197
[2671]	2198	IF (prt_level>=10) THEN
	2199	PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
	2200	wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
	2201	ENDIF
	2202
	2203
[2635]	2204	! -----------------------------------------------------------------
	2205	! Get back to tendencies per second
[1992]	2206
[2635]	2207	DO k = 1, klev
	2208	DO i = 1, klon
	2209
	2210	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
[2759]	2211	!jyg<
	2212	!! IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
	2213	IF (ok_qx_qw(i)) THEN
	2214	!>jyg
[2635]	2215	! cc
	2216	dtls(i, k) = dtls(i, k)/dtime
	2217	dqls(i, k) = dqls(i, k)/dtime
	2218	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
	2219	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
	2220	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
	2221	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
	2222	! c $ ,death_rate(i)*sigmaw(i)
	2223	END IF
	2224	END DO
	2225	END DO
	2226	!jyg<
	2227	DO i = 1, klon
	2228	d_sigmaw2(i) = d_sigmaw2(i)/dtime
[3208]	2229	d_awdens2(i) = d_awdens2(i)/dtime
[2635]	2230	d_wdens2(i) = d_wdens2(i)/dtime
	2231	ENDDO
	2232	!>jyg
	2233
	2234
	2235
[1992]	2236	RETURN
	2237	END SUBROUTINE wake
	2238
	2239	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
	2240	d_deltaqw, sigmaw, d_sigmaw, alpha)
	2241	! ------------------------------------------------------
	2242	! Dtermination du coefficient alpha tel que les tendances
	2243	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
	2244	! a une humidite positive dans la zone (x) et dans la zone (w).
	2245	! ------------------------------------------------------
[2197]	2246	IMPLICIT NONE
[1992]	2247
	2248	! Input
	2249	REAL qe(nlon, nl), d_qe(nlon, nl)
	2250	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
	2251	REAL sigmaw(nlon), d_sigmaw(nlon)
	2252	LOGICAL wk_adv(nlon)
	2253	INTEGER nl, nlon
	2254	! Output
	2255	REAL alpha(nlon)
	2256	! Internal variables
	2257	REAL zeta(nlon, nl)
	2258	REAL alpha1(nlon)
	2259	REAL x, a, b, c, discrim
	2260	REAL epsilon
	2261	! DATA epsilon/1.e-15/
[2197]	2262	INTEGER i,k
[1992]	2263
	2264	DO k = 1, nl
	2265	DO i = 1, nlon
	2266	IF (wk_adv(i)) THEN
	2267	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
	2268	zeta(i, k) = 0.
[1146]	2269	ELSE
[1992]	2270	zeta(i, k) = 1.
[1146]	2271	END IF
[1992]	2272	END IF
	2273	END DO
	2274	DO i = 1, nlon
	2275	IF (wk_adv(i)) THEN
	2276	x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
[2635]	2277	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i) &
	2278	(deltaqw(i,k)+d_deltaqw(i,k))
[1992]	2279	a = -d_sigmaw(i)*d_deltaqw(i, k)
	2280	b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
	2281	deltaqw(i, k)*d_sigmaw(i)
	2282	c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
	2283	discrim = bb - 4.a*c
	2284	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
	2285	IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
	2286	alpha1(i) = 1.
[1146]	2287	ELSE
[1992]	2288	IF (x>=0.) THEN
	2289	alpha1(i) = 1.
	2290	ELSE
	2291	IF (a>0.) THEN
[2635]	2292	alpha1(i) = 0.9min( (2.c)/(-b+sqrt(discrim)), &
	2293	(-b+sqrt(discrim))/(2.*a) )
[1992]	2294	ELSE IF (a==0.) THEN
	2295	alpha1(i) = 0.9*(-c/b)
	2296	ELSE
	2297	! print*,'a,b,c discrim',a,b,c discrim
[2635]	2298	alpha1(i) = 0.9max( (2.c)/(-b+sqrt(discrim)), &
	2299	(-b+sqrt(discrim))/(2.*a))
[1992]	2300	END IF
	2301	END IF
	2302	END IF
	2303	alpha(i) = min(alpha(i), alpha1(i))
	2304	END IF
	2305	END DO
	2306	END DO
[1146]	2307
[1992]	2308	RETURN
	2309	END SUBROUTINE wake_vec_modulation
[974]	2310
[879]	2311
[2635]	2312
	2313

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

Download in other formats:

Original Format