Context Navigation

wake.F90 @ 5440

Last change on this file since 5440 was 2787, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2727:2785 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 62.0 KB

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

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

Download in other formats:

Original Format