Context Navigation

wake.F90 @ 5225

Last change on this file since 5225 was 2925, checked in by fcheruy, 7 years ago
Update tree branch to trunk version
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.4 KB

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

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

Download in other formats:

Original Format