Context Navigation

wake.F90 @ 5507

Last change on this file since 5507 was 3899, checked in by lguez, 4 years ago
Old work on cold pool death.
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: 64.6 KB

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

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

Download in other formats:

Original Format