Context Navigation

wake.F90 @ 1992

Last change on this file since 1992 was 1992, checked in by lguez, 10 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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: 77.9 KB

Rev	Line
[1992]	1
[1403]	2	! $Id: wake.F90 1992 2014-03-05 13:19:12Z lguez $
[879]	3
[1992]	4	SUBROUTINE wake(p, ph, pi, dtime, sigd_con, te0, qe0, omgb, dtdwn, dqdwn, &
	5	amdwn, amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, deltatw, deltaqw, &
	6	dth, hw, sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, dp_omgb, &
	7	wdens, tu, qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
	8	d_deltat_gw, d_deltatw2, d_deltaqw2)
[1146]	9
[974]	10
[1992]	11	! **************************************************************
	12	! *
	13	! WAKE *
	14	! retour a un Pupper fixe *
	15	! *
	16	! written by : GRANDPEIX Jean-Yves 09/03/2000 *
	17	! modified by : ROEHRIG Romain 01/29/2007 *
	18	! **************************************************************
[974]	19
[1992]	20	USE dimphy
	21	IMPLICIT NONE
	22	! ============================================================================
[974]	23
	24
[1992]	25	! But : Decrire le comportement des poches froides apparaissant dans les
	26	! grands systemes convectifs, et fournir l'energie disponible pour
	27	! le declenchement de nouvelles colonnes convectives.
[974]	28
[1992]	29	! Variables d'etat : deltatw : ecart de temperature wake-undisturbed
	30	! area
	31	! deltaqw : ecart d'humidite wake-undisturbed area
	32	! sigmaw : fraction d'aire occupee par la poche.
[974]	33
[1992]	34	! Variable de sortie :
[974]	35
[1992]	36	! wape : WAke Potential Energy
	37	! fip : Front Incident Power (W/m2) - ALP
	38	! gfl : Gust Front Length per unit area (m-1)
	39	! dtls : large scale temperature tendency due to wake
	40	! dqls : large scale humidity tendency due to wake
	41	! hw : hauteur de la poche
	42	! dp_omgb : vertical gradient of large scale omega
	43	! wdens : densite de poches
	44	! omgbdth: flux of Delta_Theta transported by LS omega
	45	! dtKE : differential heating (wake - unpertubed)
	46	! dqKE : differential moistening (wake - unpertubed)
	47	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
	48	! dp_deltomg : vertical gradient of omg (s-1)
	49	! spread : spreading term in dt_wake and dq_wake
	50	! deltatw : updated temperature difference (T_w-T_u).
	51	! deltaqw : updated humidity difference (q_w-q_u).
	52	! sigmaw : updated wake fractional area.
	53	! d_deltat_gw : delta T tendency due to GW
[974]	54
[1992]	55	! Variables d'entree :
[974]	56
[1992]	57	! aire : aire de la maille
	58	! te0 : temperature dans l'environnement (K)
	59	! qe0 : humidite dans l'environnement (kg/kg)
	60	! omgb : vitesse verticale moyenne sur la maille (Pa/s)
	61	! dtdwn: source de chaleur due aux descentes (K/s)
	62	! dqdwn: source d'humidite due aux descentes (kg/kg/s)
	63	! dta : source de chaleur due courants satures et detrain (K/s)
	64	! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
	65	! amdwn: flux de masse total des descentes, par unite de
	66	! surface de la maille (kg/m2/s)
	67	! amup : flux de masse total des ascendances, par unite de
	68	! surface de la maille (kg/m2/s)
	69	! p : pressions aux milieux des couches (Pa)
	70	! ph : pressions aux interfaces (Pa)
	71	! pi : (p/p_0)**kapa (adim)
	72	! dtime: increment temporel (s)
[974]	73
[1992]	74	! Variables internes :
[974]	75
[1992]	76	! rhow : masse volumique de la poche froide
	77	! rho : environment density at P levels
	78	! rhoh : environment density at Ph levels
	79	! te : environment temperature \| may change within
	80	! qe : environment humidity \| sub-time-stepping
	81	! the : environment potential temperature
	82	! thu : potential temperature in undisturbed area
	83	! tu : temperature in undisturbed area
	84	! qu : humidity in undisturbed area
	85	! dp_omgb: vertical gradient og LS omega
	86	! omgbw : wake average vertical omega
	87	! dp_omgbw: vertical gradient of omgbw
	88	! omgbdq : flux of Delta_q transported by LS omega
	89	! dth : potential temperature diff. wake-undist.
	90	! th1 : first pot. temp. for vertical advection (=thu)
	91	! th2 : second pot. temp. for vertical advection (=thw)
	92	! q1 : first humidity for vertical advection
	93	! q2 : second humidity for vertical advection
	94	! d_deltatw : terme de redistribution pour deltatw
	95	! d_deltaqw : terme de redistribution pour deltaqw
	96	! deltatw0 : deltatw initial
	97	! deltaqw0 : deltaqw initial
	98	! hw0 : hw initial
	99	! sigmaw0: sigmaw initial
	100	! amflux : horizontal mass flux through wake boundary
	101	! wdens_ref: initial number of wakes per unit area (3D) or per
	102	! unit length (2D), at the beginning of each time step
	103	! Tgw : 1 sur la période de onde de gravité
	104	! Cgw : vitesse de propagation de onde de gravité
	105	! LL : distance entre 2 poches
[974]	106
[1992]	107	! -------------------------------------------------------------------------
	108	! Déclaration de variables
	109	! -------------------------------------------------------------------------
[1146]	110
[1992]	111	include "dimensions.h"
	112	include "YOMCST.h"
	113	include "cvthermo.h"
	114	include "iniprint.h"
[974]	115
[1992]	116	! Arguments en entree
	117	! --------------------
[974]	118
[1992]	119	REAL, DIMENSION (klon, klev) :: p, pi
	120	REAL, DIMENSION (klon, klev+1) :: ph, omgb
	121	REAL dtime
	122	REAL, DIMENSION (klon, klev) :: te0, qe0
	123	REAL, DIMENSION (klon, klev) :: dtdwn, dqdwn
	124	REAL, DIMENSION (klon, klev) :: wdtpbl, wdqpbl
	125	REAL, DIMENSION (klon, klev) :: udtpbl, udqpbl
	126	REAL, DIMENSION (klon, klev) :: amdwn, amup
	127	REAL, DIMENSION (klon, klev) :: dta, dqa
	128	REAL, DIMENSION (klon) :: sigd_con
[974]	129
[1992]	130	! Sorties
	131	! --------
[974]	132
[1992]	133	REAL, DIMENSION (klon, klev) :: deltatw, deltaqw, dth
	134	REAL, DIMENSION (klon, klev) :: tu, qu
	135	REAL, DIMENSION (klon, klev) :: dtls, dqls
	136	REAL, DIMENSION (klon, klev) :: dtke, dqke
	137	REAL, DIMENSION (klon, klev) :: dtpbl, dqpbl
	138	REAL, DIMENSION (klon, klev) :: spread
	139	REAL, DIMENSION (klon, klev) :: d_deltatgw
	140	REAL, DIMENSION (klon, klev) :: d_deltatw2, d_deltaqw2
	141	REAL, DIMENSION (klon, klev+1) :: omgbdth, omg
	142	REAL, DIMENSION (klon, klev) :: dp_omgb, dp_deltomg
	143	REAL, DIMENSION (klon, klev) :: d_deltat_gw
	144	REAL, DIMENSION (klon) :: hw, sigmaw, wape, fip, gfl, cstar
	145	REAL, DIMENSION (klon) :: wdens
	146	INTEGER, DIMENSION (klon) :: ktopw
[974]	147
[1992]	148	! Variables internes
	149	! -------------------
[974]	150
[1992]	151	! Variables à fixer
	152	REAL alon
	153	REAL coefgw
	154	REAL :: wdens_ref
	155	REAL stark
	156	REAL alpk
	157	REAL delta_t_min
	158	INTEGER nsub
	159	REAL dtimesub
	160	REAL sigmad, hwmin, wapecut
	161	REAL :: sigmaw_max
	162	REAL :: dens_rate
	163	REAL wdens0
	164	! IM 080208
	165	LOGICAL, DIMENSION (klon) :: gwake
[974]	166
[1992]	167	! Variables de sauvegarde
	168	REAL, DIMENSION (klon, klev) :: deltatw0
	169	REAL, DIMENSION (klon, klev) :: deltaqw0
	170	REAL, DIMENSION (klon, klev) :: te, qe
	171	REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
[974]	172
[1992]	173	! Variables pour les GW
	174	REAL, DIMENSION (klon) :: ll
	175	REAL, DIMENSION (klon, klev) :: n2
	176	REAL, DIMENSION (klon, klev) :: cgw
	177	REAL, DIMENSION (klon, klev) :: tgw
[1403]	178
[1992]	179	! Variables liées au calcul de hw
	180	REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
	181	REAL, DIMENSION (klon) :: sum_dth
	182	REAL, DIMENSION (klon) :: dthmin
	183	REAL, DIMENSION (klon) :: z, dz, hw0
	184	INTEGER, DIMENSION (klon) :: ktop, kupper
[1403]	185
[1992]	186	! Sub-timestep tendencies and related variables
	187	REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
	188	REAL d_te(klon, klev), d_qe(klon, klev)
	189	REAL d_sigmaw(klon), alpha(klon)
	190	REAL q0_min(klon), q1_min(klon)
	191	LOGICAL wk_adv(klon), ok_qx_qw(klon)
	192	REAL epsilon
	193	DATA epsilon/1.E-15/
[974]	194
[1992]	195	! Autres variables internes
	196	INTEGER isubstep, k, i
[974]	197
[1992]	198	REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
	199	REAL, DIMENSION (klon) :: sum_dq, sum_rho
	200	REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
	201	REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
	202	REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
	203	REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
[974]	204
[1992]	205	REAL, DIMENSION (klon, klev) :: rho, rhow
	206	REAL, DIMENSION (klon, klev+1) :: rhoh
	207	REAL, DIMENSION (klon, klev) :: rhow_moyen
	208	REAL, DIMENSION (klon, klev) :: zh
	209	REAL, DIMENSION (klon, klev+1) :: zhh
	210	REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
[974]	211
[1992]	212	REAL, DIMENSION (klon, klev) :: the, thu
[974]	213
[1992]	214	! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
[974]	215
[1992]	216	REAL, DIMENSION (klon, klev+1) :: omgbw
	217	REAL, DIMENSION (klon) :: pupper
	218	REAL, DIMENSION (klon) :: omgtop
	219	REAL, DIMENSION (klon, klev) :: dp_omgbw
	220	REAL, DIMENSION (klon) :: ztop, dztop
	221	REAL, DIMENSION (klon, klev) :: alpha_up
[974]	222
[1992]	223	REAL, DIMENSION (klon) :: rre1, rre2
	224	REAL :: rrd1, rrd2
	225	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
	226	REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
	227	REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
	228	REAL, DIMENSION (klon, klev) :: omgbdq
[974]	229
[1992]	230	REAL, DIMENSION (klon) :: ff, gg
	231	REAL, DIMENSION (klon) :: wape2, cstar2, heff
[974]	232
[1992]	233	REAL, DIMENSION (klon, klev) :: crep
	234	REAL crep_upper, crep_sol
[974]	235
[1992]	236	REAL, DIMENSION (klon, klev) :: ppi
[974]	237
[1992]	238	! cc nrlmd
	239	REAL, DIMENSION (klon) :: death_rate, nat_rate
	240	REAL, DIMENSION (klon, klev) :: entr
	241	REAL, DIMENSION (klon, klev) :: detr
[974]	242
[1992]	243	! -------------------------------------------------------------------------
	244	! Initialisations
	245	! -------------------------------------------------------------------------
[974]	246
[1992]	247	! print*, 'wake initialisations'
[974]	248
[1992]	249	! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
	250	! -------------------------------------------------------------------------
[974]	251
[1992]	252	DATA wapecut, sigmad, hwmin/5., .02, 10./
	253	! cc nrlmd
	254	DATA sigmaw_max/0.4/
	255	DATA dens_rate/0.1/
	256	! cc
	257	! Longueur de maille (en m)
	258	! -------------------------------------------------------------------------
[974]	259
[1992]	260	! ALON = 3.e5
	261	alon = 1.E6
[974]	262
	263
[1992]	264	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
[974]	265
[1992]	266	! coefgw : Coefficient pour les ondes de gravité
	267	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
	268	! wdens : Densité de poche froide par maille
	269	! -------------------------------------------------------------------------
[974]	270
[1992]	271	! cc nrlmd coefgw=10
	272	! coefgw=1
	273	! wdens0 = 1.0/(alon**2)
	274	! cc nrlmd wdens = 1.0/(alon**2)
	275	! cc nrlmd stark = 0.50
	276	! CRtest
	277	! cc nrlmd alpk=0.1
	278	! alpk = 1.0
	279	! alpk = 0.5
	280	! alpk = 0.05
[1146]	281
[1992]	282	stark = 0.33
	283	alpk = 0.25
	284	wdens_ref = 8.E-12
	285	coefgw = 4.
	286	crep_upper = 0.9
	287	crep_sol = 1.0
[974]	288
[1992]	289	! cc nrlmd Lecture du fichier wake_param.data
	290	OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999)
	291	READ (99, *, END=9998) stark
	292	READ (99, *, END=9998) alpk
	293	READ (99, *, END=9998) wdens_ref
	294	READ (99, *, END=9998) coefgw
	295	9998 CONTINUE
	296	CLOSE (99)
	297	9999 CONTINUE
[974]	298
[1992]	299	! Initialisation de toutes des densites a wdens_ref.
	300	! Les densites peuvent evoluer si les poches debordent
	301	! (voir au tout debut de la boucle sur les substeps)
	302	wdens = wdens_ref
[974]	303
[1992]	304	! print*,'stark',stark
	305	! print*,'alpk',alpk
	306	! print*,'wdens',wdens
	307	! print*,'coefgw',coefgw
	308	! cc
	309	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
	310	! -------------------------------------------------------------------------
[974]	311
[1992]	312	delta_t_min = 0.2
[974]	313
[1992]	314	! 1. - Save initial values and initialize tendencies
	315	! --------------------------------------------------
[974]	316
[1992]	317	DO k = 1, klev
	318	DO i = 1, klon
	319	ppi(i, k) = pi(i, k)
	320	deltatw0(i, k) = deltatw(i, k)
	321	deltaqw0(i, k) = deltaqw(i, k)
	322	te(i, k) = te0(i, k)
	323	qe(i, k) = qe0(i, k)
	324	dtls(i, k) = 0.
	325	dqls(i, k) = 0.
	326	d_deltat_gw(i, k) = 0.
	327	d_te(i, k) = 0.
	328	d_qe(i, k) = 0.
	329	d_deltatw(i, k) = 0.
	330	d_deltaqw(i, k) = 0.
	331	! IM 060508 beg
	332	d_deltatw2(i, k) = 0.
	333	d_deltaqw2(i, k) = 0.
	334	! IM 060508 end
	335	END DO
	336	END DO
	337	! sigmaw1=sigmaw
	338	! IF (sigd_con.GT.sigmaw1) THEN
	339	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
	340	! ENDIF
	341	DO i = 1, klon
	342	! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
	343	sigmaw(i) = amax1(sigmaw(i), sigmad)
	344	sigmaw(i) = amin1(sigmaw(i), 0.99)
	345	sigmaw0(i) = sigmaw(i)
	346	wape(i) = 0.
	347	wape2(i) = 0.
	348	d_sigmaw(i) = 0.
	349	ktopw(i) = 0
	350	END DO
[974]	351
	352
[1992]	353	! 2. - Prognostic part
	354	! --------------------
[974]	355
	356
[1992]	357	! 2.1 - Undisturbed area and Wake integrals
	358	! ---------------------------------------------------------
[974]	359
[1992]	360	DO i = 1, klon
	361	z(i) = 0.
	362	ktop(i) = 0
	363	kupper(i) = 0
	364	sum_thu(i) = 0.
	365	sum_tu(i) = 0.
	366	sum_qu(i) = 0.
	367	sum_thvu(i) = 0.
	368	sum_dth(i) = 0.
	369	sum_dq(i) = 0.
	370	sum_rho(i) = 0.
	371	sum_dtdwn(i) = 0.
	372	sum_dqdwn(i) = 0.
[974]	373
[1992]	374	av_thu(i) = 0.
	375	av_tu(i) = 0.
	376	av_qu(i) = 0.
	377	av_thvu(i) = 0.
	378	av_dth(i) = 0.
	379	av_dq(i) = 0.
	380	av_rho(i) = 0.
	381	av_dtdwn(i) = 0.
	382	av_dqdwn(i) = 0.
	383	END DO
[974]	384
[1992]	385	! Distance between wakes
	386	DO i = 1, klon
	387	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
	388	END DO
	389	! Potential temperatures and humidity
	390	! ----------------------------------------------------------
	391	DO k = 1, klev
	392	DO i = 1, klon
	393	! write(,)'wake 1',i,k,rd,te(i,k)
	394	rho(i, k) = p(i, k)/(rd*te(i,k))
	395	! write(,)'wake 2',rho(i,k)
	396	IF (k==1) THEN
	397	! write(,)'wake 3',i,k,rd,te(i,k)
	398	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
	399	! write(,)'wake 4',i,k,rd,te(i,k)
	400	zhh(i, k) = 0
	401	ELSE
	402	! write(,)'wake 5',rd,(te(i,k)+te(i,k-1))
	403	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
	404	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
	405	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
	406	END IF
	407	! write(,)'wake 7',ppi(i,k)
	408	the(i, k) = te(i, k)/ppi(i, k)
	409	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
	410	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
	411	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
	412	! write(,)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
	413	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
	414	dth(i, k) = deltatw(i, k)/ppi(i, k)
	415	END DO
	416	END DO
[1403]	417
[1992]	418	DO k = 1, klev - 1
	419	DO i = 1, klon
	420	IF (k==1) THEN
	421	n2(i, k) = 0
	422	ELSE
	423	n2(i, k) = amax1(0., -rg*2/the(i,k)rho(i,k)*(the(i,k+1)-the(i, &
	424	k-1))/(p(i,k+1)-p(i,k-1)))
	425	END IF
	426	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
[1403]	427
[1992]	428	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
	429	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
	430	END DO
	431	END DO
[974]	432
[1992]	433	DO i = 1, klon
	434	n2(i, klev) = 0
	435	zh(i, klev) = 0
	436	cgw(i, klev) = 0
	437	tgw(i, klev) = 0
	438	END DO
[974]	439
[1992]	440	! Calcul de la masse volumique moyenne de la colonne (bdlmd)
	441	! -----------------------------------------------------------------
[974]	442
[1992]	443	DO k = 1, klev
	444	DO i = 1, klon
	445	epaisseur1(i, k) = 0.
	446	epaisseur2(i, k) = 0.
	447	END DO
	448	END DO
[974]	449
[1992]	450	DO i = 1, klon
	451	epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
	452	epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
	453	rhow_moyen(i, 1) = rhow(i, 1)
	454	END DO
[974]	455
[1992]	456	DO k = 2, klev
	457	DO i = 1, klon
	458	epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
	459	epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
	460	rhow_moyen(i, k) = (rhow_moyen(i,k-1)epaisseur2(i,k-1)+rhow(i,k) &
	461	epaisseur1(i,k))/epaisseur2(i, k)
	462	END DO
	463	END DO
[974]	464
	465
[1992]	466	! Choose an integration bound well above wake top
	467	! -----------------------------------------------------------------
[974]	468
[1992]	469	! Pupper = 50000. ! melting level
	470	! Pupper = 60000.
	471	! Pupper = 80000. ! essais pour case_e
	472	DO i = 1, klon
	473	pupper(i) = 0.6*ph(i, 1)
	474	pupper(i) = max(pupper(i), 45000.)
	475	! cc Pupper(i) = 60000.
	476	END DO
[1146]	477
[1403]	478
[1992]	479	! Determine Wake top pressure (Ptop) from buoyancy integral
	480	! --------------------------------------------------------
[1403]	481
[1992]	482	! -1/ Pressure of the level where dth becomes less than delta_t_min.
	483
	484	DO i = 1, klon
	485	ptop_provis(i) = ph(i, 1)
	486	END DO
	487	DO k = 2, klev
	488	DO i = 1, klon
	489
	490	! IM v3JYG; ptop_provis(i).LT. ph(i,1)
	491
	492	IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
	493	ptop_provis(i)==ph(i,1)) THEN
	494	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	495	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	496	END IF
	497	END DO
	498	END DO
	499
	500	! -2/ dth integral
	501
	502	DO i = 1, klon
	503	sum_dth(i) = 0.
	504	dthmin(i) = -delta_t_min
	505	z(i) = 0.
	506	END DO
	507
	508	DO k = 1, klev
	509	DO i = 1, klon
	510	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
	511	IF (dz(i)>0) THEN
	512	z(i) = z(i) + dz(i)
	513	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	514	dthmin(i) = amin1(dthmin(i), dth(i,k))
	515	END IF
	516	END DO
	517	END DO
	518
	519	! -3/ height of triangle with area= sum_dth and base = dthmin
	520
	521	DO i = 1, klon
	522	hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
	523	hw0(i) = amax1(hwmin, hw0(i))
	524	END DO
	525
	526	! -4/ now, get Ptop
	527
	528	DO i = 1, klon
	529	z(i) = 0.
	530	ptop(i) = ph(i, 1)
	531	END DO
	532
	533	DO k = 1, klev
	534	DO i = 1, klon
	535	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
	536	IF (dz(i)>0) THEN
	537	z(i) = z(i) + dz(i)
	538	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
	539	END IF
	540	END DO
	541	END DO
	542
	543
	544	! -5/ Determination de ktop et kupper
	545
	546	DO k = klev, 1, -1
	547	DO i = 1, klon
	548	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	549	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
	550	END DO
	551	END DO
	552
	553	! On evite kupper = 1 et kupper = klev
	554	DO i = 1, klon
	555	kupper(i) = max(kupper(i), 2)
	556	kupper(i) = min(kupper(i), klev-1)
	557	END DO
	558
	559
	560	! -6/ Correct ktop and ptop
	561
	562	DO i = 1, klon
	563	ptop_new(i) = ptop(i)
	564	END DO
	565	DO k = klev, 2, -1
	566	DO i = 1, klon
	567	IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
	568	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
	569	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	570	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	571	END IF
	572	END DO
	573	END DO
	574
	575	DO i = 1, klon
	576	ptop(i) = ptop_new(i)
	577	END DO
	578
	579	DO k = klev, 1, -1
	580	DO i = 1, klon
	581	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	582	END DO
	583	END DO
	584
	585	! -5/ Set deltatw & deltaqw to 0 above kupper
	586
	587	DO k = 1, klev
	588	DO i = 1, klon
	589	IF (k>=kupper(i)) THEN
	590	deltatw(i, k) = 0.
	591	deltaqw(i, k) = 0.
	592	END IF
	593	END DO
	594	END DO
	595
	596
	597	! Vertical gradient of LS omega
	598
	599	DO k = 1, klev
	600	DO i = 1, klon
	601	IF (k<=kupper(i)) THEN
	602	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
	603	END IF
	604	END DO
	605	END DO
	606
	607	! Integrals (and wake top level number)
	608	! --------------------------------------
	609
	610	! Initialize sum_thvu to 1st level virt. pot. temp.
	611
	612	DO i = 1, klon
	613	z(i) = 1.
	614	dz(i) = 1.
	615	sum_thvu(i) = thu(i, 1)(1.+epsqu(i,1))*dz(i)
	616	sum_dth(i) = 0.
	617	END DO
	618
	619	DO k = 1, klev
	620	DO i = 1, klon
	621	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	622	IF (dz(i)>0) THEN
	623	z(i) = z(i) + dz(i)
	624	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	625	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	626	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
	627	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsqu(i,k))*dz(i)
	628	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	629	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	630	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	631	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	632	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	633	END IF
	634	END DO
	635	END DO
	636
	637	DO i = 1, klon
	638	hw0(i) = z(i)
	639	END DO
	640
	641
	642	! 2.1 - WAPE and mean forcing computation
	643	! ---------------------------------------
	644
	645	! ---------------------------------------
	646
	647	! Means
	648
	649	DO i = 1, klon
	650	av_thu(i) = sum_thu(i)/hw0(i)
	651	av_tu(i) = sum_tu(i)/hw0(i)
	652	av_qu(i) = sum_qu(i)/hw0(i)
	653	av_thvu(i) = sum_thvu(i)/hw0(i)
	654	! av_thve = sum_thve/hw0
	655	av_dth(i) = sum_dth(i)/hw0(i)
	656	av_dq(i) = sum_dq(i)/hw0(i)
	657	av_rho(i) = sum_rho(i)/hw0(i)
	658	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	659	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	660
	661	wape(i) = -rghw0(i)(av_dth(i)+eps(av_thu(i)av_dq(i)+av_dth(i)*av_qu(i &
	662	)+av_dth(i)*av_dq(i)))/av_thvu(i)
	663	END DO
	664
	665	! 2.2 Prognostic variable update
	666	! ------------------------------
	667
	668	! Filter out bad wakes
	669
	670	DO k = 1, klev
	671	DO i = 1, klon
	672	IF (wape(i)<0.) THEN
	673	deltatw(i, k) = 0.
	674	deltaqw(i, k) = 0.
	675	dth(i, k) = 0.
	676	END IF
	677	END DO
	678	END DO
	679
	680	DO i = 1, klon
	681	IF (wape(i)<0.) THEN
	682	wape(i) = 0.
	683	cstar(i) = 0.
	684	hw(i) = hwmin
	685	sigmaw(i) = amax1(sigmad, sigd_con(i))
	686	fip(i) = 0.
	687	gwake(i) = .FALSE.
	688	ELSE
	689	cstar(i) = starksqrt(2.wape(i))
	690	gwake(i) = .TRUE.
	691	END IF
	692	END DO
	693
	694
	695	! Check qx and qw positivity
	696	! --------------------------
	697	DO i = 1, klon
	698	q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, &
	699	1)+(1.-sigmaw(i))*deltaqw(i,1)))
	700	END DO
	701	DO k = 2, klev
	702	DO i = 1, klon
	703	q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, &
	704	k)+(1.-sigmaw(i))*deltaqw(i,k)))
	705	IF (q1_min(i)<=q0_min(i)) THEN
	706	q0_min(i) = q1_min(i)
	707	END IF
	708	END DO
	709	END DO
	710
	711	DO i = 1, klon
	712	ok_qx_qw(i) = q0_min(i) >= 0.
	713	alpha(i) = 1.
	714	END DO
	715
	716	! C -----------------------------------------------------------------
	717	! Sub-time-stepping
	718	! -----------------
	719
	720	nsub = 10
	721	dtimesub = dtime/nsub
	722
	723	! ------------------------------------------------------------
	724	DO isubstep = 1, nsub
	725	! ------------------------------------------------------------
	726
	727	! wk_adv is the logical flag enabling wake evolution in the time advance
	728	! loop
	729	DO i = 1, klon
	730	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
	731	END DO
	732
	733	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
	734	! négatif de ktop à kupper --------
	735	! cc On calcule pour cela une densité wdens0 pour laquelle on
	736	! aurait un entrainement nul ---
	737	DO i = 1, klon
	738	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
	739	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
	740	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
	741	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
	742	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
	743	wdens0 = (sigmaw(i)/(4.3.14))((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( &
	744	ph(i,1)-pupper(i))cstar(i)))*(2)
	745	IF (wdens(i)<=wdens0*1.1) THEN
	746	wdens(i) = wdens0
	747	END IF
	748	! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
	749	! c $ ,ph(i,1)-pupper(i)',
	750	! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
	751	! c $ ,ph(i,1)-pupper(i)
	752	END IF
	753	END DO
	754
	755	! cc nrlmd
	756
	757	DO i = 1, klon
	758	IF (wk_adv(i)) THEN
[1403]	759	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
[1992]	760	sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
	761	END IF
	762	END DO
	763	DO i = 1, klon
	764	IF (wk_adv(i)) THEN
	765	! cc nrlmd Introduction du taux de mortalité des poches et
	766	! test sur sigmaw_max=0.4
	767	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
	768	IF (sigmaw(i)>=sigmaw_max) THEN
	769	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
	770	ELSE
	771	death_rate(i) = 0.
	772	END IF
	773	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
	774	dtimesub
	775	! $ - nat_rate(i)sigmaw(i)dtimesub
	776	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
	777	! c $ death_rate(i),ktop(i),kupper(i)',
	778	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
	779	! c $ death_rate(i),ktop(i),kupper(i)
[1403]	780
[1992]	781	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
	782	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
	783	! wdens = wdens0/(10.*sigmaw)
	784	! sigmaw =max(sigmaw,sigd_con)
	785	! sigmaw =max(sigmaw,sigmad)
	786	END IF
	787	END DO
	788
	789
	790	! calcul de la difference de vitesse verticale poche - zone non perturbee
	791	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
	792	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on
	793	! definit
	794	! IM 060208 au niveau k=1..?
	795	DO k = 1, klev
	796	DO i = 1, klon
	797	IF (wk_adv(i)) THEN !!! nrlmd
	798	dp_deltomg(i, k) = 0.
	799	END IF
	800	END DO
	801	END DO
	802	DO k = 1, klev + 1
	803	DO i = 1, klon
	804	IF (wk_adv(i)) THEN !!! nrlmd
	805	omg(i, k) = 0.
	806	END IF
	807	END DO
	808	END DO
	809
	810	DO i = 1, klon
	811	IF (wk_adv(i)) THEN
	812	z(i) = 0.
	813	omg(i, 1) = 0.
	814	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
	815	END IF
	816	END DO
	817
	818	DO k = 2, klev
	819	DO i = 1, klon
	820	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
[974]	821	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
[1992]	822	z(i) = z(i) + dz(i)
	823	dp_deltomg(i, k) = dp_deltomg(i, 1)
	824	omg(i, k) = dp_deltomg(i, 1)*z(i)
	825	END IF
	826	END DO
	827	END DO
	828
	829	DO i = 1, klon
	830	IF (wk_adv(i)) THEN
	831	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
	832	ztop(i) = z(i) + dztop(i)
	833	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
	834	END IF
	835	END DO
	836
	837	! -----------------
	838	! From m/s to Pa/s
	839	! -----------------
	840
	841	DO i = 1, klon
	842	IF (wk_adv(i)) THEN
	843	omgtop(i) = -rho(i, ktop(i))rgomgtop(i)
	844	dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
	845	END IF
	846	END DO
	847
	848	DO k = 1, klev
	849	DO i = 1, klon
	850	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
	851	omg(i, k) = -rho(i, k)rgomg(i, k)
	852	dp_deltomg(i, k) = dp_deltomg(i, 1)
	853	END IF
	854	END DO
	855	END DO
	856
	857	! raccordement lineaire de omg de ptop a pupper
	858
	859	DO i = 1, klon
	860	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
	861	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
	862	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
	863	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
	864	(ptop(i)-pupper(i))
	865	END IF
	866	END DO
	867
	868	! c DO i=1,klon
	869	! c print*,'Pente entre 0 et kupper (référence)'
	870	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
	871	! c print*,'Pente entre ktop et kupper'
	872	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
	873	! c ENDDO
	874	! c
	875	DO k = 1, klev
	876	DO i = 1, klon
	877	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
	878	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
	879	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
	880	END IF
	881	END DO
	882	END DO
	883	! cc nrlmd
	884	! c DO i=1,klon
	885	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
	886	! c END DO
	887	! cc
	888
	889
	890	! -- Compute wake average vertical velocity omgbw
	891
	892
	893	DO k = 1, klev + 1
	894	DO i = 1, klon
[1146]	895	IF (wk_adv(i)) THEN
[1992]	896	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
	897	END IF
	898	END DO
	899	END DO
	900	! -- and its vertical gradient dp_omgbw
	901
	902	DO k = 1, klev
	903	DO i = 1, klon
[1146]	904	IF (wk_adv(i)) THEN
[1992]	905	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
	906	END IF
	907	END DO
	908	END DO
[974]	909
[1992]	910	! -- Upstream coefficients for omgb velocity
	911	! -- (alpha_up(k) is the coefficient of the value at level k)
	912	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
	913	DO k = 1, klev
	914	DO i = 1, klon
	915	IF (wk_adv(i)) THEN
	916	alpha_up(i, k) = 0.
	917	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
	918	END IF
	919	END DO
	920	END DO
[974]	921
[1992]	922	! Matrix expressing [The,deltatw] from [Th1,Th2]
[974]	923
[1992]	924	DO i = 1, klon
	925	IF (wk_adv(i)) THEN
	926	rre1(i) = 1. - sigmaw(i)
	927	rre2(i) = sigmaw(i)
	928	END IF
	929	END DO
	930	rrd1 = -1.
	931	rrd2 = 1.
[974]	932
[1992]	933	! -- Get [Th1,Th2], dth and [q1,q2]
[974]	934
[1992]	935	DO k = 1, klev
	936	DO i = 1, klon
	937	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
	938	dth(i, k) = deltatw(i, k)/ppi(i, k)
	939	th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
	940	th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
	941	q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
	942	q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
	943	END IF
	944	END DO
	945	END DO
[974]	946
[1992]	947	DO i = 1, klon
	948	IF (wk_adv(i)) THEN !!! nrlmd
	949	d_th1(i, 1) = 0.
	950	d_th2(i, 1) = 0.
	951	d_dth(i, 1) = 0.
	952	d_q1(i, 1) = 0.
	953	d_q2(i, 1) = 0.
	954	d_dq(i, 1) = 0.
	955	END IF
	956	END DO
[974]	957
[1992]	958	DO k = 2, klev
	959	DO i = 1, klon
	960	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
	961	d_th1(i, k) = th1(i, k-1) - th1(i, k)
	962	d_th2(i, k) = th2(i, k-1) - th2(i, k)
	963	d_dth(i, k) = dth(i, k-1) - dth(i, k)
	964	d_q1(i, k) = q1(i, k-1) - q1(i, k)
	965	d_q2(i, k) = q2(i, k-1) - q2(i, k)
	966	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
	967	END IF
	968	END DO
	969	END DO
[1146]	970
[1992]	971	DO i = 1, klon
	972	IF (wk_adv(i)) THEN
	973	omgbdth(i, 1) = 0.
	974	omgbdq(i, 1) = 0.
	975	END IF
	976	END DO
[1277]	977
[1992]	978	DO k = 2, klev
	979	DO i = 1, klon
	980	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
	981	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
	982	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
	983	END IF
	984	END DO
	985	END DO
[1403]	986
[1992]	987	! -----------------------------------------------------------------
	988	DO k = 1, klev
	989	DO i = 1, klon
	990	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
	991	! -----------------------------------------------------------------
[974]	992
[1992]	993	! Compute redistribution (advective) term
[1403]	994
[1992]	995	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
	996	(rrd1omg(i,k)sigmaw(i)d_th1(i,k)-rrd2omg(i,k+1)*(1.-sigmaw( &
	997	i))d_th2(i,k+1)-(1.-alpha_up(i,k))omgbdth(i,k)-alpha_up(i,k+1)* &
	998	omgbdth(i,k+1))*ppi(i, k)
	999	! print*,'d_deltatw=',d_deltatw(i,k)
[1403]	1000
[1992]	1001	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
	1002	(rrd1omg(i,k)sigmaw(i)d_q1(i,k)-rrd2omg(i,k+1)*(1.-sigmaw( &
	1003	i))d_q2(i,k+1)-(1.-alpha_up(i,k))omgbdq(i,k)-alpha_up(i,k+1)* &
	1004	omgbdq(i,k+1))
	1005	! print*,'d_deltaqw=',d_deltaqw(i,k)
[974]	1006
[1992]	1007	! and increment large scale tendencies
[974]	1008
	1009
	1010
	1011
[1992]	1012	! C
	1013	! -----------------------------------------------------------------
	1014	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i, &
	1015	k)-rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
	1016	k+1)) & ! cc nrlmd $
	1017	! -sigmaw(i)(1.-sigmaw(i))dth(i,k)*dp_deltomg(i,k)
	1018	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
	1019	k)-ph(i,k+1)) & ! cc
	1020	)*ppi(i, k)
[974]	1021
[1992]	1022	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i, &
	1023	k)-rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
	1024	k+1)) & ! cc nrlmd $
	1025	! -sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*dp_deltomg(i,k)
	1026	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i, &
	1027	k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
	1028	)
	1029	! cc nrlmd
	1030	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
	1031	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i, &
	1032	k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k)
[1403]	1033
[1992]	1034	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i, &
	1035	k)/(ph(i,k)-ph(i,k+1))))
[1403]	1036
[1992]	1037	END IF
	1038	! cc
	1039	END DO
	1040	END DO
	1041	! ------------------------------------------------------------------
[974]	1042
[1992]	1043	! Increment state variables
[974]	1044
[1992]	1045	DO k = 1, klev
	1046	DO i = 1, klon
	1047	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
	1048	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1049	! cc
[974]	1050
[1146]	1051
[974]	1052
[1992]	1053	! Coefficient de répartition
[974]	1054
[1992]	1055	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
	1056	(ph(i,kupper(i))-ph(i,1))
	1057	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i &
	1058	,kupper(i)))
[974]	1059
	1060
[1992]	1061	! Reintroduce compensating subsidence term.
[1146]	1062
[1992]	1063	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
	1064	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
	1065	! . /(1-sigmaw)
	1066	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
	1067	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
	1068	! . /(1-sigmaw)
[974]	1069
[1992]	1070	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
	1071	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
	1072	! . /(1-sigmaw)
	1073	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
	1074	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
	1075	! . /(1-sigmaw)
[974]	1076
[1992]	1077	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
	1078	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
	1079	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
[974]	1080
[1992]	1081	dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
	1082	dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
	1083	! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
[1146]	1084
[1992]	1085	! cc nrlmd Prise en compte du taux de mortalité
	1086	! cc Définitions de entr, detr
	1087	detr(i, k) = 0.
[1146]	1088
[1992]	1089	entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
	1090	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
[1146]	1091
[1992]	1092	spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
	1093	! cc spread(i,k) =
	1094	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
	1095	! cc $ sigmaw(i)
[1146]	1096
	1097
[1992]	1098	! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
	1099	! Jingmei
[1146]	1100
[1992]	1101	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
	1102	! & Tgw(i,k),deltatw(i,k)
	1103	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
	1104	dtimesub
	1105	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
	1106	ff(i) = d_deltatw(i, k)/dtimesub
[1403]	1107
[1992]	1108	! Sans GW
[1403]	1109
[1992]	1110	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
[974]	1111
[1992]	1112	! GW formule 1
	1113
	1114	! deltatw(k) = deltatw(k)+dtimesub*
	1115	! $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
	1116
	1117	! GW formule 2
	1118
	1119	IF (dtimesub*tgw(i,k)<1.E-10) THEN
	1120	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
	1121	! $
	1122	! -spread(i,k)*deltatw(i,k)
	1123	-entr(i,k)deltatw(i,k)/sigmaw(i)-(death_rate(i)sigmaw( &
	1124	i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
	1125	-tgw(i,k)*deltatw(i,k))
	1126	ELSE
	1127	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i, &
	1128	k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $
	1129	! -spread(i,k)*deltatw(i,k)
	1130	-entr(i,k)deltatw(i,k)/sigmaw(i)-(death_rate(i)sigmaw( &
	1131	i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
	1132	-tgw(i,k)*deltatw(i,k))
	1133	END IF
	1134
	1135	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1136
	1137	gg(i) = d_deltaqw(i, k)/dtimesub
	1138
	1139	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $
	1140	! -spread(i,k)*deltaqw(i,k))
	1141	-entr(i,k)deltaqw(i,k)/sigmaw(i)-(death_rate(i)sigmaw(i)+detr( &
	1142	i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
	1143	! cc
	1144
	1145	! cc nrlmd
	1146	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
	1147	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
	1148	! cc
	1149	END IF
	1150	END DO
	1151	END DO
	1152
	1153
	1154	! Scale tendencies so that water vapour remains positive in w and x.
	1155
	1156	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
	1157	d_deltaqw, sigmaw, d_sigmaw, alpha)
	1158
	1159	! cc nrlmd
	1160	! c print*,'alpha'
	1161	! c do i=1,klon
	1162	! c print*,alpha(i)
	1163	! c end do
	1164	! cc
	1165	DO k = 1, klev
	1166	DO i = 1, klon
	1167	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1168	d_te(i, k) = alpha(i)*d_te(i, k)
	1169	d_qe(i, k) = alpha(i)*d_qe(i, k)
	1170	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
	1171	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
	1172	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
	1173	END IF
	1174	END DO
	1175	END DO
	1176	DO i = 1, klon
	1177	IF (wk_adv(i)) THEN
	1178	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
	1179	END IF
	1180	END DO
	1181
	1182	! Update large scale variables and wake variables
	1183	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
	1184	! IM 060208 DO k = 1,kupper(i)
	1185	DO k = 1, klev
	1186	DO i = 1, klon
	1187	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1188	dtls(i, k) = dtls(i, k) + d_te(i, k)
	1189	dqls(i, k) = dqls(i, k) + d_qe(i, k)
	1190	! cc nrlmd
	1191	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
	1192	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
	1193	! cc
	1194	END IF
	1195	END DO
	1196	END DO
	1197	DO k = 1, klev
	1198	DO i = 1, klon
	1199	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1200	te(i, k) = te0(i, k) + dtls(i, k)
	1201	qe(i, k) = qe0(i, k) + dqls(i, k)
	1202	the(i, k) = te(i, k)/ppi(i, k)
	1203	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
	1204	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
	1205	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1206	! c print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
	1207	! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
	1208	END IF
	1209	END DO
	1210	END DO
	1211	DO i = 1, klon
	1212	IF (wk_adv(i)) THEN
	1213	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
	1214	END IF
	1215	END DO
	1216
	1217
	1218	! Determine Ptop from buoyancy integral
	1219	! ---------------------------------------
	1220
	1221	! - 1/ Pressure of the level where dth changes sign.
	1222
	1223	DO i = 1, klon
	1224	IF (wk_adv(i)) THEN
	1225	ptop_provis(i) = ph(i, 1)
	1226	END IF
	1227	END DO
	1228
	1229	DO k = 2, klev
	1230	DO i = 1, klon
	1231	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
	1232	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
	1233	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	1234	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	1235	END IF
	1236	END DO
	1237	END DO
	1238
	1239	! - 2/ dth integral
	1240
	1241	DO i = 1, klon
	1242	IF (wk_adv(i)) THEN !!! nrlmd
	1243	sum_dth(i) = 0.
	1244	dthmin(i) = -delta_t_min
[974]	1245	z(i) = 0.
[1992]	1246	END IF
	1247	END DO
	1248
	1249	DO k = 1, klev
	1250	DO i = 1, klon
	1251	IF (wk_adv(i)) THEN
	1252	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
	1253	IF (dz(i)>0) THEN
	1254	z(i) = z(i) + dz(i)
	1255	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1256	dthmin(i) = amin1(dthmin(i), dth(i,k))
	1257	END IF
	1258	END IF
	1259	END DO
	1260	END DO
	1261
	1262	! - 3/ height of triangle with area= sum_dth and base = dthmin
	1263
	1264	DO i = 1, klon
	1265	IF (wk_adv(i)) THEN
	1266	hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
	1267	hw(i) = amax1(hwmin, hw(i))
	1268	END IF
	1269	END DO
	1270
	1271	! - 4/ now, get Ptop
	1272
	1273	DO i = 1, klon
	1274	IF (wk_adv(i)) THEN !!! nrlmd
	1275	ktop(i) = 0
	1276	z(i) = 0.
	1277	END IF
	1278	END DO
	1279
	1280	DO k = 1, klev
	1281	DO i = 1, klon
	1282	IF (wk_adv(i)) THEN
	1283	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
	1284	IF (dz(i)>0) THEN
	1285	z(i) = z(i) + dz(i)
	1286	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
	1287	ktop(i) = k
	1288	END IF
	1289	END IF
	1290	END DO
	1291	END DO
	1292
	1293	! 4.5/Correct ktop and ptop
	1294
	1295	DO i = 1, klon
	1296	IF (wk_adv(i)) THEN
	1297	ptop_new(i) = ptop(i)
	1298	END IF
	1299	END DO
	1300
	1301	DO k = klev, 2, -1
	1302	DO i = 1, klon
	1303	! IM v3JYG; IF (k .GE. ktop(i)
	1304	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
	1305	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
	1306	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	1307	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	1308	END IF
	1309	END DO
	1310	END DO
	1311
	1312
	1313	DO i = 1, klon
	1314	IF (wk_adv(i)) THEN
	1315	ptop(i) = ptop_new(i)
	1316	END IF
	1317	END DO
	1318
	1319	DO k = klev, 1, -1
	1320	DO i = 1, klon
	1321	IF (wk_adv(i)) THEN !!! nrlmd
	1322	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	1323	END IF
	1324	END DO
	1325	END DO
	1326
	1327	! 5/ Set deltatw & deltaqw to 0 above kupper
	1328
	1329	DO k = 1, klev
	1330	DO i = 1, klon
	1331	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
	1332	deltatw(i, k) = 0.
	1333	deltaqw(i, k) = 0.
	1334	END IF
	1335	END DO
	1336	END DO
	1337
	1338
	1339	! -------------Cstar computation---------------------------------
	1340	DO i = 1, klon
	1341	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1342	sum_thu(i) = 0.
	1343	sum_tu(i) = 0.
	1344	sum_qu(i) = 0.
	1345	sum_thvu(i) = 0.
	1346	sum_dth(i) = 0.
	1347	sum_dq(i) = 0.
	1348	sum_rho(i) = 0.
	1349	sum_dtdwn(i) = 0.
	1350	sum_dqdwn(i) = 0.
	1351
	1352	av_thu(i) = 0.
[1992]	1353	av_tu(i) = 0.
	1354	av_qu(i) = 0.
[974]	1355	av_thvu(i) = 0.
	1356	av_dth(i) = 0.
	1357	av_dq(i) = 0.
[1992]	1358	av_rho(i) = 0.
	1359	av_dtdwn(i) = 0.
[974]	1360	av_dqdwn(i) = 0.
[1992]	1361	END IF
	1362	END DO
[974]	1363
[1992]	1364	! Integrals (and wake top level number)
	1365	! --------------------------------------
[974]	1366
[1992]	1367	! Initialize sum_thvu to 1st level virt. pot. temp.
[974]	1368
[1992]	1369	DO i = 1, klon
	1370	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1371	z(i) = 1.
	1372	dz(i) = 1.
[1992]	1373	sum_thvu(i) = thu(i, 1)(1.+epsqu(i,1))*dz(i)
[974]	1374	sum_dth(i) = 0.
[1992]	1375	END IF
	1376	END DO
[974]	1377
[1992]	1378	DO k = 1, klev
	1379	DO i = 1, klon
	1380	IF (wk_adv(i)) THEN !!! nrlmd
	1381	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	1382	IF (dz(i)>0) THEN
	1383	z(i) = z(i) + dz(i)
	1384	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	1385	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	1386	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
	1387	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsqu(i,k))*dz(i)
	1388	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1389	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	1390	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	1391	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	1392	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	1393	END IF
	1394	END IF
	1395	END DO
	1396	END DO
	1397
	1398	DO i = 1, klon
	1399	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1400	hw0(i) = z(i)
[1992]	1401	END IF
	1402	END DO
[974]	1403
	1404
[1992]	1405	! - WAPE and mean forcing computation
	1406	! ---------------------------------------
	1407
	1408	! ---------------------------------------
	1409
	1410	! Means
	1411
	1412	DO i = 1, klon
	1413	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1414	av_thu(i) = sum_thu(i)/hw0(i)
	1415	av_tu(i) = sum_tu(i)/hw0(i)
	1416	av_qu(i) = sum_qu(i)/hw0(i)
	1417	av_thvu(i) = sum_thvu(i)/hw0(i)
	1418	av_dth(i) = sum_dth(i)/hw0(i)
	1419	av_dq(i) = sum_dq(i)/hw0(i)
	1420	av_rho(i) = sum_rho(i)/hw0(i)
	1421	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	1422	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	1423
[1992]	1424	wape(i) = -rghw0(i)(av_dth(i)+eps(av_thu(i)av_dq(i)+av_dth(i)* &
	1425	av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
	1426	END IF
	1427	END DO
[974]	1428
[1992]	1429	! Filter out bad wakes
[974]	1430
[1992]	1431	DO k = 1, klev
	1432	DO i = 1, klon
	1433	IF (wk_adv(i)) THEN !!! nrlmd
	1434	IF (wape(i)<0.) THEN
	1435	deltatw(i, k) = 0.
	1436	deltaqw(i, k) = 0.
	1437	dth(i, k) = 0.
	1438	END IF
	1439	END IF
	1440	END DO
	1441	END DO
[974]	1442
[1992]	1443	DO i = 1, klon
	1444	IF (wk_adv(i)) THEN !!! nrlmd
	1445	IF (wape(i)<0.) THEN
	1446	wape(i) = 0.
	1447	cstar(i) = 0.
	1448	hw(i) = hwmin
	1449	sigmaw(i) = max(sigmad, sigd_con(i))
	1450	fip(i) = 0.
	1451	gwake(i) = .FALSE.
	1452	ELSE
	1453	cstar(i) = starksqrt(2.wape(i))
	1454	gwake(i) = .TRUE.
	1455	END IF
	1456	END IF
	1457	END DO
	1458
	1459	END DO ! end sub-timestep loop
	1460
	1461	! -----------------------------------------------------------------
	1462	! Get back to tendencies per second
	1463
	1464	DO k = 1, klev
	1465	DO i = 1, klon
	1466
	1467	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
	1468	IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
	1469	! cc
	1470	dtls(i, k) = dtls(i, k)/dtime
	1471	dqls(i, k) = dqls(i, k)/dtime
	1472	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
	1473	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
	1474	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
	1475	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
	1476	! c $ ,death_rate(i)*sigmaw(i)
	1477	END IF
	1478	END DO
	1479	END DO
	1480
	1481
	1482	! ----------------------------------------------------------
	1483	! Determine wake final state; recompute wape, cstar, ktop;
	1484	! filter out bad wakes.
	1485	! ----------------------------------------------------------
	1486
	1487	! 2.1 - Undisturbed area and Wake integrals
	1488	! ---------------------------------------------------------
	1489
	1490	DO i = 1, klon
	1491	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
	1492	IF (ok_qx_qw(i)) THEN
	1493	! cc
	1494	z(i) = 0.
	1495	sum_thu(i) = 0.
	1496	sum_tu(i) = 0.
	1497	sum_qu(i) = 0.
	1498	sum_thvu(i) = 0.
	1499	sum_dth(i) = 0.
	1500	sum_dq(i) = 0.
	1501	sum_rho(i) = 0.
	1502	sum_dtdwn(i) = 0.
	1503	sum_dqdwn(i) = 0.
	1504
	1505	av_thu(i) = 0.
	1506	av_tu(i) = 0.
	1507	av_qu(i) = 0.
	1508	av_thvu(i) = 0.
	1509	av_dth(i) = 0.
	1510	av_dq(i) = 0.
	1511	av_rho(i) = 0.
	1512	av_dtdwn(i) = 0.
	1513	av_dqdwn(i) = 0.
	1514	END IF
	1515	END DO
	1516	! Potential temperatures and humidity
	1517	! ----------------------------------------------------------
	1518
	1519	DO k = 1, klev
	1520	DO i = 1, klon
	1521	! cc nrlmd IF ( wk_adv(i)) THEN
	1522	IF (ok_qx_qw(i)) THEN
	1523	! cc
	1524	rho(i, k) = p(i, k)/(rd*te(i,k))
	1525	IF (k==1) THEN
	1526	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
	1527	zhh(i, k) = 0
	1528	ELSE
	1529	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
	1530	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
	1531	END IF
	1532	the(i, k) = te(i, k)/ppi(i, k)
	1533	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
	1534	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
	1535	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
	1536	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
	1537	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1538	END IF
	1539	END DO
	1540	END DO
	1541
	1542	! Integrals (and wake top level number)
	1543	! -----------------------------------------------------------
	1544
	1545	! Initialize sum_thvu to 1st level virt. pot. temp.
	1546
	1547	DO i = 1, klon
	1548	! cc nrlmd IF ( wk_adv(i)) THEN
	1549	IF (ok_qx_qw(i)) THEN
	1550	! cc
	1551	z(i) = 1.
	1552	dz(i) = 1.
	1553	sum_thvu(i) = thu(i, 1)(1.+epsqu(i,1))*dz(i)
	1554	sum_dth(i) = 0.
	1555	END IF
	1556	END DO
	1557
	1558	DO k = 1, klev
	1559	DO i = 1, klon
	1560	! cc nrlmd IF ( wk_adv(i)) THEN
	1561	IF (ok_qx_qw(i)) THEN
	1562	! cc
	1563	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	1564	IF (dz(i)>0) THEN
	1565	z(i) = z(i) + dz(i)
	1566	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	1567	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	1568	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
	1569	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsqu(i,k))*dz(i)
	1570	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1571	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	1572	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	1573	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	1574	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	1575	END IF
	1576	END IF
	1577	END DO
	1578	END DO
	1579
	1580	DO i = 1, klon
	1581	! cc nrlmd IF ( wk_adv(i)) THEN
	1582	IF (ok_qx_qw(i)) THEN
	1583	! cc
	1584	hw0(i) = z(i)
	1585	END IF
	1586	END DO
	1587
	1588	! - WAPE and mean forcing computation
	1589	! -------------------------------------------------------------
	1590
	1591	! Means
	1592
	1593	DO i = 1, klon
	1594	! cc nrlmd IF ( wk_adv(i)) THEN
	1595	IF (ok_qx_qw(i)) THEN
	1596	! cc
	1597	av_thu(i) = sum_thu(i)/hw0(i)
	1598	av_tu(i) = sum_tu(i)/hw0(i)
	1599	av_qu(i) = sum_qu(i)/hw0(i)
	1600	av_thvu(i) = sum_thvu(i)/hw0(i)
	1601	av_dth(i) = sum_dth(i)/hw0(i)
	1602	av_dq(i) = sum_dq(i)/hw0(i)
	1603	av_rho(i) = sum_rho(i)/hw0(i)
	1604	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	1605	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	1606
	1607	wape2(i) = -rghw0(i)(av_dth(i)+eps(av_thu(i)av_dq(i)+av_dth(i)* &
	1608	av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
	1609	END IF
	1610	END DO
	1611
	1612	! Prognostic variable update
	1613	! ------------------------------------------------------------
	1614
	1615	! Filter out bad wakes
	1616
	1617	DO k = 1, klev
	1618	DO i = 1, klon
	1619	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
	1620	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
	1621	! cc
	1622	deltatw(i, k) = 0.
	1623	deltaqw(i, k) = 0.
	1624	dth(i, k) = 0.
	1625	END IF
	1626	END DO
	1627	END DO
	1628
	1629
	1630	DO i = 1, klon
	1631	! cc nrlmd IF ( wk_adv(i)) THEN
	1632	IF (ok_qx_qw(i)) THEN
	1633	! cc
	1634	IF (wape2(i)<0.) THEN
[974]	1635	wape2(i) = 0.
[1992]	1636	cstar2(i) = 0.
[974]	1637	hw(i) = hwmin
[1992]	1638	sigmaw(i) = amax1(sigmad, sigd_con(i))
[974]	1639	fip(i) = 0.
	1640	gwake(i) = .FALSE.
	1641	ELSE
[1992]	1642	IF (prt_level>=10) PRINT *, 'wape2>0'
	1643	cstar2(i) = starksqrt(2.wape2(i))
[974]	1644	gwake(i) = .TRUE.
[1992]	1645	END IF
	1646	END IF
	1647	END DO
[974]	1648
[1992]	1649	DO i = 1, klon
	1650	! cc nrlmd IF ( wk_adv(i)) THEN
	1651	IF (ok_qx_qw(i)) THEN
	1652	! cc
	1653	ktopw(i) = ktop(i)
	1654	END IF
	1655	END DO
[974]	1656
[1992]	1657	DO i = 1, klon
	1658	! cc nrlmd IF ( wk_adv(i)) THEN
	1659	IF (ok_qx_qw(i)) THEN
	1660	! cc
	1661	IF (ktopw(i)>0 .AND. gwake(i)) THEN
[1403]	1662
[1992]	1663	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
	1664	! cc heff = 600.
	1665	! Utilisation de la hauteur hw
	1666	! c heff = 0.7*hw
	1667	heff(i) = hw(i)
[1403]	1668
[1992]	1669	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
	1670	sqrt(sigmaw(i)wdens(i)3.14)
	1671	fip(i) = alpk*fip(i)
	1672	! jyg2
	1673	ELSE
	1674	fip(i) = 0.
	1675	END IF
	1676	END IF
	1677	END DO
[1146]	1678
[1992]	1679	! Limitation de sigmaw
	1680
	1681	! cc nrlmd
	1682	! DO i=1,klon
	1683	! IF (OK_qx_qw(i)) THEN
	1684	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
	1685	! ENDIF
	1686	! ENDDO
	1687	! cc
	1688	DO k = 1, klev
	1689	DO i = 1, klon
	1690
	1691	! cc nrlmd On maintient désormais constant sigmaw en régime
	1692	! permanent
	1693	! cc IF ((sigmaw(i).GT.sigmaw_max).or.
	1694	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
	1695	.NOT. ok_qx_qw(i)) THEN
	1696	! cc
	1697	dtls(i, k) = 0.
	1698	dqls(i, k) = 0.
	1699	deltatw(i, k) = 0.
	1700	deltaqw(i, k) = 0.
	1701	END IF
	1702	END DO
	1703	END DO
	1704
	1705	! cc nrlmd On maintient désormais constant sigmaw en régime permanent
	1706	DO i = 1, klon
	1707	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
	1708	.NOT. ok_qx_qw(i)) THEN
	1709	wape(i) = 0.
	1710	cstar(i) = 0.
	1711	hw(i) = hwmin
	1712	sigmaw(i) = sigmad
	1713	fip(i) = 0.
	1714	ELSE
	1715	wape(i) = wape2(i)
	1716	cstar(i) = cstar2(i)
	1717	END IF
	1718	! c print*,'wape wape2 ktopw OK_qx_qw =',
	1719	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
	1720	END DO
	1721
	1722
	1723	RETURN
	1724	END SUBROUTINE wake
	1725
	1726	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
	1727	d_deltaqw, sigmaw, d_sigmaw, alpha)
	1728	! ------------------------------------------------------
	1729	! Dtermination du coefficient alpha tel que les tendances
	1730	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
	1731	! a une humidite positive dans la zone (x) et dans la zone (w).
	1732	! ------------------------------------------------------
	1733
	1734
	1735	! Input
	1736	REAL qe(nlon, nl), d_qe(nlon, nl)
	1737	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
	1738	REAL sigmaw(nlon), d_sigmaw(nlon)
	1739	LOGICAL wk_adv(nlon)
	1740	INTEGER nl, nlon
	1741	! Output
	1742	REAL alpha(nlon)
	1743	! Internal variables
	1744	REAL zeta(nlon, nl)
	1745	REAL alpha1(nlon)
	1746	REAL x, a, b, c, discrim
	1747	REAL epsilon
	1748	! DATA epsilon/1.e-15/
	1749
	1750	DO k = 1, nl
	1751	DO i = 1, nlon
	1752	IF (wk_adv(i)) THEN
	1753	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
	1754	zeta(i, k) = 0.
[1146]	1755	ELSE
[1992]	1756	zeta(i, k) = 1.
[1146]	1757	END IF
[1992]	1758	END IF
	1759	END DO
	1760	DO i = 1, nlon
	1761	IF (wk_adv(i)) THEN
	1762	x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
	1763	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i)(deltaqw(i,k)+ &
	1764	d_deltaqw(i,k))
	1765	a = -d_sigmaw(i)*d_deltaqw(i, k)
	1766	b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
	1767	deltaqw(i, k)*d_sigmaw(i)
	1768	c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
	1769	discrim = bb - 4.a*c
	1770	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
	1771	IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
	1772	alpha1(i) = 1.
[1146]	1773	ELSE
[1992]	1774	IF (x>=0.) THEN
	1775	alpha1(i) = 1.
	1776	ELSE
	1777	IF (a>0.) THEN
	1778	alpha1(i) = 0.9min((2.c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
	1779	))/(2.*a))
	1780	ELSE IF (a==0.) THEN
	1781	alpha1(i) = 0.9*(-c/b)
	1782	ELSE
	1783	! print*,'a,b,c discrim',a,b,c discrim
	1784	alpha1(i) = 0.9max((2.c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
	1785	))/(2.*a))
	1786	END IF
	1787	END IF
	1788	END IF
	1789	alpha(i) = min(alpha(i), alpha1(i))
	1790	END IF
	1791	END DO
	1792	END DO
[1146]	1793
[1992]	1794	RETURN
	1795	END SUBROUTINE wake_vec_modulation
[974]	1796
[1992]	1797	SUBROUTINE wake_scal(p, ph, ppi, dtime, sigd_con, te0, qe0, omgb, dtdwn, &
	1798	dqdwn, amdwn, amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, deltatw, &
	1799	deltaqw, dth, hw, sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, &
	1800	dp_omgb, wdens, tu, qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, &
	1801	spread, cstar, d_deltat_gw, d_deltatw2, d_deltaqw2)
[879]	1802
[1992]	1803	! **************************************************************
	1804	! *
	1805	! WAKE *
	1806	! retour a un Pupper fixe *
	1807	! *
	1808	! written by : GRANDPEIX Jean-Yves 09/03/2000 *
	1809	! modified by : ROEHRIG Romain 01/29/2007 *
	1810	! **************************************************************
[879]	1811
[1992]	1812	USE dimphy
	1813	IMPLICIT NONE
	1814	! ============================================================================
[879]	1815
	1816
[1992]	1817	! But : Decrire le comportement des poches froides apparaissant dans les
	1818	! grands systemes convectifs, et fournir l'energie disponible pour
	1819	! le declenchement de nouvelles colonnes convectives.
[879]	1820
[1992]	1821	! Variables d'etat : deltatw : ecart de temperature wake-undisturbed
	1822	! area
	1823	! deltaqw : ecart d'humidite wake-undisturbed area
	1824	! sigmaw : fraction d'aire occupee par la poche.
[879]	1825
[1992]	1826	! Variable de sortie :
[879]	1827
[1992]	1828	! wape : WAke Potential Energy
	1829	! fip : Front Incident Power (W/m2) - ALP
	1830	! gfl : Gust Front Length per unit area (m-1)
	1831	! dtls : large scale temperature tendency due to wake
	1832	! dqls : large scale humidity tendency due to wake
	1833	! hw : hauteur de la poche
	1834	! dp_omgb : vertical gradient of large scale omega
	1835	! omgbdth: flux of Delta_Theta transported by LS omega
	1836	! dtKE : differential heating (wake - unpertubed)
	1837	! dqKE : differential moistening (wake - unpertubed)
	1838	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
	1839	! dp_deltomg : vertical gradient of omg (s-1)
	1840	! spread : spreading term in dt_wake and dq_wake
	1841	! deltatw : updated temperature difference (T_w-T_u).
	1842	! deltaqw : updated humidity difference (q_w-q_u).
	1843	! sigmaw : updated wake fractional area.
	1844	! d_deltat_gw : delta T tendency due to GW
[879]	1845
[1992]	1846	! Variables d'entree :
[879]	1847
[1992]	1848	! aire : aire de la maille
	1849	! te0 : temperature dans l'environnement (K)
	1850	! qe0 : humidite dans l'environnement (kg/kg)
	1851	! omgb : vitesse verticale moyenne sur la maille (Pa/s)
	1852	! dtdwn: source de chaleur due aux descentes (K/s)
	1853	! dqdwn: source d'humidite due aux descentes (kg/kg/s)
	1854	! dta : source de chaleur due courants satures et detrain (K/s)
	1855	! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
	1856	! amdwn: flux de masse total des descentes, par unite de
	1857	! surface de la maille (kg/m2/s)
	1858	! amup : flux de masse total des ascendances, par unite de
	1859	! surface de la maille (kg/m2/s)
	1860	! p : pressions aux milieux des couches (Pa)
	1861	! ph : pressions aux interfaces (Pa)
	1862	! ppi : (p/p_0)**kapa (adim)
	1863	! dtime: increment temporel (s)
[879]	1864
[1992]	1865	! Variables internes :
[879]	1866
[1992]	1867	! rhow : masse volumique de la poche froide
	1868	! rho : environment density at P levels
	1869	! rhoh : environment density at Ph levels
	1870	! te : environment temperature \| may change within
	1871	! qe : environment humidity \| sub-time-stepping
	1872	! the : environment potential temperature
	1873	! thu : potential temperature in undisturbed area
	1874	! tu : temperature in undisturbed area
	1875	! qu : humidity in undisturbed area
	1876	! dp_omgb: vertical gradient og LS omega
	1877	! omgbw : wake average vertical omega
	1878	! dp_omgbw: vertical gradient of omgbw
	1879	! omgbdq : flux of Delta_q transported by LS omega
	1880	! dth : potential temperature diff. wake-undist.
	1881	! th1 : first pot. temp. for vertical advection (=thu)
	1882	! th2 : second pot. temp. for vertical advection (=thw)
	1883	! q1 : first humidity for vertical advection
	1884	! q2 : second humidity for vertical advection
	1885	! d_deltatw : terme de redistribution pour deltatw
	1886	! d_deltaqw : terme de redistribution pour deltaqw
	1887	! deltatw0 : deltatw initial
	1888	! deltaqw0 : deltaqw initial
	1889	! hw0 : hw initial
	1890	! sigmaw0: sigmaw initial
	1891	! amflux : horizontal mass flux through wake boundary
	1892	! wdens : number of wakes per unit area (3D) or per
	1893	! unit length (2D)
	1894	! Tgw : 1 sur la période de onde de gravité
	1895	! Cgw : vitesse de propagation de onde de gravité
	1896	! LL : distance entre 2 poches
[879]	1897
[1992]	1898	! -------------------------------------------------------------------------
	1899	! Déclaration de variables
	1900	! -------------------------------------------------------------------------
[879]	1901
[1992]	1902	include "dimensions.h"
	1903	! ccc include "dimphy.h"
	1904	include "YOMCST.h"
	1905	include "cvthermo.h"
	1906	include "iniprint.h"
[879]	1907
[1992]	1908	! Arguments en entree
	1909	! --------------------
[879]	1910
[1992]	1911	REAL p(klev), ph(klev+1), ppi(klev)
	1912	REAL dtime
	1913	REAL te0(klev), qe0(klev)
	1914	REAL omgb(klev+1)
	1915	REAL dtdwn(klev), dqdwn(klev)
	1916	REAL wdtpbl(klev), wdqpbl(klev)
	1917	REAL udtpbl(klev), udqpbl(klev)
	1918	REAL amdwn(klev), amup(klev)
	1919	REAL dta(klev), dqa(klev)
	1920	REAL sigd_con
[879]	1921
[1992]	1922	! Sorties
	1923	! --------
[879]	1924
[1992]	1925	REAL deltatw(klev), deltaqw(klev), dth(klev)
	1926	REAL tu(klev), qu(klev)
	1927	REAL dtls(klev), dqls(klev)
	1928	REAL dtke(klev), dqke(klev)
	1929	REAL dtpbl(klev), dqpbl(klev)
	1930	REAL spread(klev)
	1931	REAL d_deltatgw(klev)
	1932	REAL d_deltatw2(klev), d_deltaqw2(klev)
	1933	REAL omgbdth(klev+1), omg(klev+1)
	1934	REAL dp_omgb(klev), dp_deltomg(klev)
	1935	REAL d_deltat_gw(klev)
	1936	REAL hw, sigmaw, wape, fip, gfl, cstar
	1937	INTEGER ktopw
[879]	1938
[1992]	1939	! Variables internes
	1940	! -------------------
[879]	1941
[1992]	1942	! Variables à fixer
	1943	REAL alon
	1944	REAL coefgw
	1945	REAL wdens0, wdens
	1946	REAL stark
	1947	REAL alpk
	1948	REAL delta_t_min
	1949	REAL pupper
	1950	INTEGER nsub
	1951	REAL dtimesub
	1952	REAL sigmad, hwmin
[879]	1953
[1992]	1954	! Variables de sauvegarde
	1955	REAL deltatw0(klev)
	1956	REAL deltaqw0(klev)
	1957	REAL te(klev), qe(klev)
	1958	REAL sigmaw0, sigmaw1
[879]	1959
[1992]	1960	! Variables pour les GW
	1961	REAL ll
	1962	REAL n2(klev)
	1963	REAL cgw(klev)
	1964	REAL tgw(klev)
[879]	1965
[1992]	1966	! Variables liées au calcul de hw
	1967	REAL ptop_provis, ptop, ptop_new
	1968	REAL sum_dth
	1969	REAL dthmin
	1970	REAL z, dz, hw0
	1971	INTEGER ktop, kupper
[879]	1972
[1992]	1973	! Autres variables internes
	1974	INTEGER isubstep, k
[879]	1975
[1992]	1976	REAL sum_thu, sum_tu, sum_qu, sum_thvu
	1977	REAL sum_dq, sum_rho
	1978	REAL sum_dtdwn, sum_dqdwn
	1979	REAL av_thu, av_tu, av_qu, av_thvu
	1980	REAL av_dth, av_dq, av_rho
	1981	REAL av_dtdwn, av_dqdwn
[879]	1982
[1992]	1983	REAL rho(klev), rhoh(klev+1), rhow(klev)
	1984	REAL rhow_moyen(klev)
	1985	REAL zh(klev), zhh(klev+1)
	1986	REAL epaisseur1(klev), epaisseur2(klev)
[879]	1987
[1992]	1988	REAL the(klev), thu(klev)
[879]	1989
[1992]	1990	REAL d_deltatw(klev), d_deltaqw(klev)
[879]	1991
[1992]	1992	REAL omgbw(klev+1), omgtop
	1993	REAL dp_omgbw(klev)
	1994	REAL ztop, dztop
	1995	REAL alpha_up(klev)
[879]	1996
[1992]	1997	REAL rre1, rre2, rrd1, rrd2
	1998	REAL th1(klev), th2(klev), q1(klev), q2(klev)
	1999	REAL d_th1(klev), d_th2(klev), d_dth(klev)
	2000	REAL d_q1(klev), d_q2(klev), d_dq(klev)
	2001	REAL omgbdq(klev)
[879]	2002
[1992]	2003	REAL ff, gg
	2004	REAL wape2, cstar2, heff
[879]	2005
[1992]	2006	REAL crep(klev)
	2007	REAL crep_upper, crep_sol
[879]	2008
[1992]	2009	! -------------------------------------------------------------------------
	2010	! Initialisations
	2011	! -------------------------------------------------------------------------
[879]	2012
[1992]	2013	! print*, 'wake initialisations'
[879]	2014
[1992]	2015	! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
	2016	! -------------------------------------------------------------------------
[879]	2017
[1992]	2018	DATA sigmad, hwmin/.02, 10./
[879]	2019
[1992]	2020	! Longueur de maille (en m)
	2021	! -------------------------------------------------------------------------
[879]	2022
[1992]	2023	! ALON = 3.e5
	2024	alon = 1.E6
[879]	2025
	2026
[1992]	2027	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
[879]	2028
[1992]	2029	! coefgw : Coefficient pour les ondes de gravité
	2030	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
	2031	! wdens : Densité de poche froide par maille
	2032	! -------------------------------------------------------------------------
[879]	2033
[1992]	2034	coefgw = 10
	2035	! coefgw=1
	2036	! wdens0 = 1.0/(alon**2)
	2037	wdens = 1.0/(alon**2)
	2038	stark = 0.50
	2039	! CRtest
	2040	alpk = 0.1
	2041	! alpk = 1.0
	2042	! alpk = 0.5
	2043	! alpk = 0.05
	2044	crep_upper = 0.9
	2045	crep_sol = 1.0
[879]	2046
	2047
[1992]	2048	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
	2049	! -------------------------------------------------------------------------
[879]	2050
[1992]	2051	delta_t_min = 0.2
[879]	2052
	2053
[1992]	2054	! 1. - Save initial values and initialize tendencies
	2055	! --------------------------------------------------
[879]	2056
[1992]	2057	DO k = 1, klev
	2058	deltatw0(k) = deltatw(k)
	2059	deltaqw0(k) = deltaqw(k)
	2060	te(k) = te0(k)
	2061	qe(k) = qe0(k)
	2062	dtls(k) = 0.
	2063	dqls(k) = 0.
	2064	d_deltat_gw(k) = 0.
	2065	d_deltatw2(k) = 0.
	2066	d_deltaqw2(k) = 0.
	2067	END DO
	2068	! sigmaw1=sigmaw
	2069	! IF (sigd_con.GT.sigmaw1) THEN
	2070	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
	2071	! ENDIF
	2072	sigmaw = max(sigmaw, sigd_con)
	2073	sigmaw = max(sigmaw, sigmad)
	2074	sigmaw = min(sigmaw, 0.99)
	2075	sigmaw0 = sigmaw
	2076	! wdens=wdens0/(10.*sigmaw)
	2077	! IF (sigd_con.GT.sigmaw1) THEN
	2078	! print*, 'sigmaw1,sigd1', sigmaw, sigd_con
	2079	! ENDIF
[879]	2080
[1992]	2081	! 2. - Prognostic part
	2082	! =========================================================
[879]	2083
[1992]	2084	! print *, 'prognostic wake computation'
[879]	2085
	2086
[1992]	2087	! 2.1 - Undisturbed area and Wake integrals
	2088	! ---------------------------------------------------------
[879]	2089
[1992]	2090	z = 0.
	2091	ktop = 0
	2092	kupper = 0
	2093	sum_thu = 0.
	2094	sum_tu = 0.
	2095	sum_qu = 0.
	2096	sum_thvu = 0.
	2097	sum_dth = 0.
	2098	sum_dq = 0.
	2099	sum_rho = 0.
	2100	sum_dtdwn = 0.
	2101	sum_dqdwn = 0.
[879]	2102
[1992]	2103	av_thu = 0.
	2104	av_tu = 0.
	2105	av_qu = 0.
	2106	av_thvu = 0.
	2107	av_dth = 0.
	2108	av_dq = 0.
	2109	av_rho = 0.
	2110	av_dtdwn = 0.
	2111	av_dqdwn = 0.
[879]	2112
[1992]	2113	! Potential temperatures and humidity
	2114	! ----------------------------------------------------------
[879]	2115
[1992]	2116	DO k = 1, klev
	2117	rho(k) = p(k)/(rd*te(k))
	2118	IF (k==1) THEN
	2119	rhoh(k) = ph(k)/(rd*te(k))
	2120	zhh(k) = 0
	2121	ELSE
	2122	rhoh(k) = ph(k)2./(rd(te(k)+te(k-1)))
	2123	zhh(k) = (ph(k)-ph(k-1))/(-rhoh(k)*rg) + zhh(k-1)
	2124	END IF
	2125	the(k) = te(k)/ppi(k)
	2126	thu(k) = (te(k)-deltatw(k)*sigmaw)/ppi(k)
	2127	tu(k) = te(k) - deltatw(k)*sigmaw
	2128	qu(k) = qe(k) - deltaqw(k)*sigmaw
	2129	rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
	2130	dth(k) = deltatw(k)/ppi(k)
	2131	ll = (1-sqrt(sigmaw))/sqrt(wdens)
	2132	END DO
[879]	2133
[1992]	2134	DO k = 1, klev - 1
	2135	IF (k==1) THEN
	2136	n2(k) = 0
	2137	ELSE
	2138	n2(k) = max(0., -rg*2/the(k)rho(k)*(the(k+1)-the(k-1))/(p(k+ &
	2139	1)-p(k-1)))
	2140	END IF
	2141	zh(k) = (zhh(k)+zhh(k+1))/2
[879]	2142
[1992]	2143	cgw(k) = sqrt(n2(k))*zh(k)
	2144	tgw(k) = coefgw*cgw(k)/ll
	2145	END DO
[879]	2146
[1992]	2147	n2(klev) = 0
	2148	zh(klev) = 0
	2149	cgw(klev) = 0
	2150	tgw(klev) = 0
[879]	2151
[1992]	2152	! Calcul de la masse volumique moyenne de la colonne
	2153	! -----------------------------------------------------------------
[879]	2154
[1992]	2155	DO k = 1, klev
	2156	epaisseur1(k) = 0.
	2157	epaisseur2(k) = 0.
	2158	END DO
[879]	2159
[1992]	2160	epaisseur1(1) = -(ph(2)-ph(1))/(rho(1)*rg) + 1.
	2161	epaisseur2(1) = -(ph(2)-ph(1))/(rho(1)*rg) + 1.
	2162	rhow_moyen(1) = rhow(1)
[879]	2163
[1992]	2164	DO k = 2, klev
	2165	epaisseur1(k) = -(ph(k+1)-ph(k))/(rho(k)*rg) + 1.
	2166	epaisseur2(k) = epaisseur2(k-1) + epaisseur1(k)
	2167	rhow_moyen(k) = (rhow_moyen(k-1)epaisseur2(k-1)+rhow(k)epaisseur1(k))/ &
	2168	epaisseur2(k)
	2169	END DO
[879]	2170
	2171
[1992]	2172	! Choose an integration bound well above wake top
	2173	! -----------------------------------------------------------------
[879]	2174
[1992]	2175	! Pupper = 50000. ! melting level
	2176	pupper = 60000.
	2177	! Pupper = 70000.
[879]	2178
	2179
[1992]	2180	! Determine Wake top pressure (Ptop) from buoyancy integral
	2181	! -----------------------------------------------------------------
[879]	2182
[1992]	2183	! -1/ Pressure of the level where dth becomes less than delta_t_min.
[879]	2184
[1992]	2185	ptop_provis = ph(1)
	2186	DO k = 2, klev
	2187	IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
	2188	ptop_provis = ((dth(k)+delta_t_min)*p(k-1)-(dth(k- &
	2189	1)+delta_t_min)*p(k))/(dth(k)-dth(k-1))
	2190	GO TO 25
	2191	END IF
	2192	END DO
	2193	25 CONTINUE
[879]	2194
[1992]	2195	! -2/ dth integral
[879]	2196
[1992]	2197	sum_dth = 0.
	2198	dthmin = -delta_t_min
	2199	z = 0.
[879]	2200
[1992]	2201	DO k = 1, klev
	2202	dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg)
	2203	IF (dz<=0) GO TO 40
	2204	z = z + dz
	2205	sum_dth = sum_dth + dth(k)*dz
	2206	dthmin = min(dthmin, dth(k))
	2207	END DO
	2208	40 CONTINUE
[879]	2209
[1992]	2210	! -3/ height of triangle with area= sum_dth and base = dthmin
[879]	2211
[1992]	2212	hw0 = 2.*sum_dth/min(dthmin, -0.5)
	2213	hw0 = max(hwmin, hw0)
[879]	2214
[1992]	2215	! -4/ now, get Ptop
[879]	2216
[1992]	2217	z = 0.
	2218	ptop = ph(1)
[879]	2219
[1992]	2220	DO k = 1, klev
	2221	dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg), hw0-z)
	2222	IF (dz<=0) GO TO 45
	2223	z = z + dz
	2224	ptop = ph(k) - rho(k)rgdz
	2225	END DO
	2226	45 CONTINUE
[879]	2227
	2228
[1992]	2229	! -5/ Determination de ktop et kupper
[879]	2230
[1992]	2231	DO k = klev, 1, -1
	2232	IF (ph(k+1)<ptop) ktop = k
	2233	IF (ph(k+1)<pupper) kupper = k
	2234	END DO
[879]	2235
[1992]	2236	! -6/ Correct ktop and ptop
[879]	2237
[1992]	2238	ptop_new = ptop
	2239	DO k = ktop, 2, -1
	2240	IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
	2241	ptop_new = ((dth(k)+delta_t_min)p(k-1)-(dth(k-1)+delta_t_min)p(k))/ &
	2242	(dth(k)-dth(k-1))
	2243	GO TO 225
	2244	END IF
	2245	END DO
	2246	225 CONTINUE
[879]	2247
[1992]	2248	ptop = ptop_new
[879]	2249
[1992]	2250	DO k = klev, 1, -1
	2251	IF (ph(k+1)<ptop) ktop = k
	2252	END DO
[879]	2253
[1992]	2254	! Set deltatw & deltaqw to 0 above kupper
	2255	! -----------------------------------------------------------
[879]	2256
[1992]	2257	DO k = kupper, klev
	2258	deltatw(k) = 0.
	2259	deltaqw(k) = 0.
	2260	END DO
[879]	2261
	2262
[1992]	2263	! Vertical gradient of LS omega
	2264	! ------------------------------------------------------------
[879]	2265
[1992]	2266	DO k = 1, kupper
	2267	dp_omgb(k) = (omgb(k+1)-omgb(k))/(ph(k+1)-ph(k))
	2268	END DO
[879]	2269
	2270
[1992]	2271	! Integrals (and wake top level number)
	2272	! -----------------------------------------------------------
[879]	2273
[1992]	2274	! Initialize sum_thvu to 1st level virt. pot. temp.
[879]	2275
[1992]	2276	z = 1.
	2277	dz = 1.
	2278	sum_thvu = thu(1)(1.+epsqu(1))*dz
	2279	sum_dth = 0.
[879]	2280
[1992]	2281	DO k = 1, klev
	2282	dz = -(max(ph(k+1),ptop)-ph(k))/(rho(k)*rg)
	2283	IF (dz<=0) GO TO 50
	2284	z = z + dz
	2285	sum_thu = sum_thu + thu(k)*dz
	2286	sum_tu = sum_tu + tu(k)*dz
	2287	sum_qu = sum_qu + qu(k)*dz
	2288	sum_thvu = sum_thvu + thu(k)(1.+epsqu(k))*dz
	2289	sum_dth = sum_dth + dth(k)*dz
	2290	sum_dq = sum_dq + deltaqw(k)*dz
	2291	sum_rho = sum_rho + rhow(k)*dz
	2292	sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
	2293	sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
	2294	END DO
	2295	50 CONTINUE
[879]	2296
[1992]	2297	hw0 = z
[879]	2298
[1992]	2299	! 2.1 - WAPE and mean forcing computation
	2300	! -------------------------------------------------------------
[879]	2301
[1992]	2302	! Means
[879]	2303
[1992]	2304	av_thu = sum_thu/hw0
	2305	av_tu = sum_tu/hw0
	2306	av_qu = sum_qu/hw0
	2307	av_thvu = sum_thvu/hw0
	2308	! av_thve = sum_thve/hw0
	2309	av_dth = sum_dth/hw0
	2310	av_dq = sum_dq/hw0
	2311	av_rho = sum_rho/hw0
	2312	av_dtdwn = sum_dtdwn/hw0
	2313	av_dqdwn = sum_dqdwn/hw0
[879]	2314
[1992]	2315	wape = -rghw0(av_dth+eps(av_thuav_dq+av_dthav_qu+av_dthav_dq))/ &
	2316	av_thvu
[879]	2317
[1992]	2318	! 2.2 Prognostic variable update
	2319	! ------------------------------------------------------------
[879]	2320
[1992]	2321	! Filter out bad wakes
[879]	2322
[1992]	2323	IF (wape<0.) THEN
	2324	IF (prt_level>=10) PRINT *, 'wape<0'
	2325	wape = 0.
	2326	hw = hwmin
	2327	sigmaw = max(sigmad, sigd_con)
	2328	fip = 0.
	2329	DO k = 1, klev
	2330	deltatw(k) = 0.
	2331	deltaqw(k) = 0.
	2332	dth(k) = 0.
	2333	END DO
	2334	ELSE
	2335	IF (prt_level>=10) PRINT *, 'wape>0'
	2336	cstar = starksqrt(2.wape)
	2337	END IF
[879]	2338
[1992]	2339	! ------------------------------------------------------------------
	2340	! Sub-time-stepping
	2341	! ------------------------------------------------------------------
[879]	2342
[1992]	2343	! nsub=36
	2344	nsub = 10
	2345	dtimesub = dtime/nsub
[879]	2346
[1992]	2347	! ------------------------------------------------------------
	2348	DO isubstep = 1, nsub
	2349	! ------------------------------------------------------------
[879]	2350
[1992]	2351	! print*,'---------------','substep=',isubstep,'-------------'
[879]	2352
[1992]	2353	! Evolution of sigmaw
[879]	2354
	2355
[1992]	2356	gfl = 2.sqrt(3.14wdens*sigmaw)
[879]	2357
[1992]	2358	sigmaw = sigmaw + gflcstardtimesub
	2359	sigmaw = min(sigmaw, 0.99) !!!!!!!!
	2360	! wdens = wdens0/(10.*sigmaw)
	2361	! sigmaw =max(sigmaw,sigd_con)
	2362	! sigmaw =max(sigmaw,sigmad)
[879]	2363
[1992]	2364	! calcul de la difference de vitesse verticale poche - zone non perturbee
[879]	2365
[1992]	2366	z = 0.
	2367	dp_deltomg(1:klev) = 0.
	2368	omg(1:klev+1) = 0.
[879]	2369
[1992]	2370	omg(1) = 0.
	2371	dp_deltomg(1) = -(gflcstar)/(sigmaw(1-sigmaw))
[879]	2372
[1992]	2373	DO k = 2, ktop
	2374	dz = -(ph(k)-ph(k-1))/(rho(k-1)*rg)
	2375	z = z + dz
	2376	dp_deltomg(k) = dp_deltomg(1)
	2377	omg(k) = dp_deltomg(1)*z
	2378	END DO
[879]	2379
[1992]	2380	dztop = -(ptop-ph(ktop))/(rho(ktop)*rg)
	2381	ztop = z + dztop
	2382	omgtop = dp_deltomg(1)*ztop
[879]	2383
	2384
[1992]	2385	! Conversion de la vitesse verticale de m/s a Pa/s
[879]	2386
[1992]	2387	omgtop = -rho(ktop)rgomgtop
	2388	dp_deltomg(1) = omgtop/(ptop-ph(1))
[879]	2389
[1992]	2390	DO k = 1, ktop
	2391	omg(k) = -rho(k)rgomg(k)
	2392	dp_deltomg(k) = dp_deltomg(1)
	2393	END DO
[879]	2394
[1992]	2395	! raccordement lineaire de omg de ptop a pupper
[879]	2396
[1992]	2397	IF (kupper>ktop) THEN
	2398	omg(kupper+1) = -rgamdwn(kupper+1)/sigmaw + rgamup(kupper+1)/(1.- &
	2399	sigmaw)
	2400	dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(ptop-pupper)
	2401	DO k = ktop + 1, kupper
	2402	dp_deltomg(k) = dp_deltomg(kupper)
	2403	omg(k) = omgtop + (ph(k)-ptop)*dp_deltomg(kupper)
	2404	END DO
	2405	END IF
[879]	2406
[1992]	2407	! Compute wake average vertical velocity omgbw
[879]	2408
[1992]	2409	DO k = 1, klev + 1
	2410	omgbw(k) = omgb(k) + (1.-sigmaw)*omg(k)
	2411	END DO
[879]	2412
[1992]	2413	! and its vertical gradient dp_omgbw
[879]	2414
[1992]	2415	DO k = 1, klev
	2416	dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
	2417	END DO
[879]	2418
	2419
[1992]	2420	! Upstream coefficients for omgb velocity
	2421	! -- (alpha_up(k) is the coefficient of the value at level k)
	2422	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
[879]	2423
[1992]	2424	DO k = 1, klev
	2425	alpha_up(k) = 0.
	2426	IF (omgb(k)>0.) alpha_up(k) = 1.
	2427	END DO
[879]	2428
[1992]	2429	! Matrix expressing [The,deltatw] from [Th1,Th2]
[879]	2430
[1992]	2431	rre1 = 1. - sigmaw
	2432	rre2 = sigmaw
	2433	rrd1 = -1.
	2434	rrd2 = 1.
[879]	2435
[1992]	2436	! Get [Th1,Th2], dth and [q1,q2]
[879]	2437
[1992]	2438	DO k = 1, kupper + 1
	2439	dth(k) = deltatw(k)/ppi(k)
	2440	th1(k) = the(k) - sigmaw*dth(k) ! undisturbed area
	2441	th2(k) = the(k) + (1.-sigmaw)*dth(k) ! wake
	2442	q1(k) = qe(k) - sigmaw*deltaqw(k) ! undisturbed area
	2443	q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
	2444	END DO
[879]	2445
[1992]	2446	d_th1(1) = 0.
	2447	d_th2(1) = 0.
	2448	d_dth(1) = 0.
	2449	d_q1(1) = 0.
	2450	d_q2(1) = 0.
	2451	d_dq(1) = 0.
[879]	2452
[1992]	2453	DO k = 2, kupper + 1 ! loop on interfaces
	2454	d_th1(k) = th1(k-1) - th1(k)
	2455	d_th2(k) = th2(k-1) - th2(k)
	2456	d_dth(k) = dth(k-1) - dth(k)
	2457	d_q1(k) = q1(k-1) - q1(k)
	2458	d_q2(k) = q2(k-1) - q2(k)
	2459	d_dq(k) = deltaqw(k-1) - deltaqw(k)
	2460	END DO
[879]	2461
[1992]	2462	omgbdth(1) = 0.
	2463	omgbdq(1) = 0.
[879]	2464
[1992]	2465	DO k = 2, kupper + 1 ! loop on interfaces
	2466	omgbdth(k) = omgb(k)*(dth(k-1)-dth(k))
	2467	omgbdq(k) = omgb(k)*(deltaqw(k-1)-deltaqw(k))
	2468	END DO
[879]	2469
	2470
[1992]	2471	! -----------------------------------------------------------------
	2472	DO k = 1, kupper - 1
	2473	! -----------------------------------------------------------------
[879]	2474
[1992]	2475	! Compute redistribution (advective) term
[879]	2476
[1992]	2477	d_deltatw(k) = dtimesub/(ph(k)-ph(k+1))(rrd1omg(k)sigmawd_th1(k)- &
	2478	rrd2omg(k+1)(1.-sigmaw)*d_th2(k+1)-(1.-alpha_up( &
	2479	k))omgbdth(k)-alpha_up(k+1)omgbdth(k+1))*ppi(k)
	2480	! print*,'d_deltatw=',d_deltatw(k)
[879]	2481
[1992]	2482	d_deltaqw(k) = dtimesub/(ph(k)-ph(k+1))(rrd1omg(k)sigmawd_q1(k)- &
	2483	rrd2omg(k+1)(1.-sigmaw)*d_q2(k+1)-(1.-alpha_up( &
	2484	k))omgbdq(k)-alpha_up(k+1)omgbdq(k+1))
	2485	! print*,'d_deltaqw=',d_deltaqw(k)
[879]	2486
[1992]	2487	! and increment large scale tendencies
[879]	2488
[1992]	2489	dtls(k) = dtls(k) + dtimesub((rre1omg(k)sigmawd_th1(k)-rre2*omg(k+ &
	2490	1)(1.-sigmaw)d_th2(k+1))/(ph(k)-ph(k+1))-sigmaw(1.-sigmaw)dth(k)* &
	2491	dp_deltomg(k))*ppi(k)
	2492	! print*,'dtls=',dtls(k)
[879]	2493
[1992]	2494	dqls(k) = dqls(k) + dtimesub((rre1omg(k)sigmawd_q1(k)-rre2*omg(k+ &
	2495	1)(1.-sigmaw)d_q2(k+1))/(ph(k)-ph(k+1))-sigmaw(1.-sigmaw)deltaqw( &
	2496	k)*dp_deltomg(k))
	2497	! print*,'dqls=',dqls(k)
[879]	2498
[1992]	2499	! -------------------------------------------------------------------
	2500	END DO
	2501	! ------------------------------------------------------------------
	2502
	2503	! Increment state variables
	2504
	2505	DO k = 1, kupper - 1
	2506
	2507	! Coefficient de répartition
	2508
	2509	crep(k) = crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
	2510	crep(k) = crep(k) + crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
	2511
	2512
	2513	! Reintroduce compensating subsidence term.
	2514
	2515	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
	2516	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
	2517	! . /(1-sigmaw)
	2518	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
	2519	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
	2520	! . /(1-sigmaw)
	2521
	2522	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
	2523	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
	2524	! . /(1-sigmaw)
	2525	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
	2526	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
	2527	! . /(1-sigmaw)
	2528
	2529	dtke(k) = (dtdwn(k)/sigmaw-dta(k)/(1.-sigmaw))
	2530	dqke(k) = (dqdwn(k)/sigmaw-dqa(k)/(1.-sigmaw))
	2531	! print*,'dtKE=',dtKE(k)
	2532	! print*,'dqKE=',dqKE(k)
	2533
	2534	dtpbl(k) = (wdtpbl(k)/sigmaw-udtpbl(k)/(1.-sigmaw))
	2535	dqpbl(k) = (wdqpbl(k)/sigmaw-udqpbl(k)/(1.-sigmaw))
	2536
	2537	spread(k) = (1.-sigmaw)dp_deltomg(k) + gflcstar/sigmaw
	2538	! print*,'spread=',spread(k)
	2539
	2540
	2541	! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
	2542	! Jingmei
	2543
	2544	d_deltat_gw(k) = d_deltat_gw(k) - tgw(k)deltatw(k)dtimesub
	2545	! print*,'d_delta_gw=',d_deltat_gw(k)
	2546	ff = d_deltatw(k)/dtimesub
	2547
	2548	! Sans GW
	2549
	2550	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
	2551
	2552	! GW formule 1
	2553
	2554	! deltatw(k) = deltatw(k)+dtimesub*
	2555	! $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
	2556
	2557	! GW formule 2
	2558
	2559	IF (dtimesub*tgw(k)<1.E-10) THEN
	2560	deltatw(k) = deltatw(k) + dtimesub(ff+dtke(k)+dtpbl(k)-spread(k) &
	2561	deltatw(k)-tgw(k)*deltatw(k))
	2562	ELSE
	2563	deltatw(k) = deltatw(k) + 1/tgw(k)(1-exp(-dtimesubtgw(k)))*(ff+dtke &
	2564	(k)+dtpbl(k)-spread(k)deltatw(k)-tgw(k)deltatw(k))
	2565	END IF
	2566
	2567	dth(k) = deltatw(k)/ppi(k)
	2568
	2569	gg = d_deltaqw(k)/dtimesub
	2570
	2571	deltaqw(k) = deltaqw(k) + dtimesub(gg+dqke(k)+dqpbl(k)-spread(k) &
	2572	deltaqw(k))
	2573
	2574	d_deltatw2(k) = d_deltatw2(k) + d_deltatw(k)
	2575	d_deltaqw2(k) = d_deltaqw2(k) + d_deltaqw(k)
	2576	END DO
	2577
	2578	! And update large scale variables
	2579
	2580	DO k = 1, kupper
	2581	te(k) = te0(k) + dtls(k)
	2582	qe(k) = qe0(k) + dqls(k)
	2583	the(k) = te(k)/ppi(k)
	2584	END DO
	2585
	2586	! Determine Ptop from buoyancy integral
	2587	! ----------------------------------------------------------------------
	2588
	2589	! -1/ Pressure of the level where dth changes sign.
	2590
	2591	ptop_provis = ph(1)
	2592
	2593	DO k = 2, klev
	2594	IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
	2595	ptop_provis = ((dth(k)+delta_t_min)*p(k-1)-(dth(k- &
	2596	1)+delta_t_min)*p(k))/(dth(k)-dth(k-1))
[879]	2597	GO TO 65
[1992]	2598	END IF
	2599	END DO
	2600	65 CONTINUE
[879]	2601
[1992]	2602	! -2/ dth integral
[879]	2603
[1992]	2604	sum_dth = 0.
	2605	dthmin = -delta_t_min
	2606	z = 0.
[879]	2607
[1992]	2608	DO k = 1, klev
	2609	dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg)
	2610	IF (dz<=0) GO TO 70
	2611	z = z + dz
	2612	sum_dth = sum_dth + dth(k)*dz
	2613	dthmin = min(dthmin, dth(k))
	2614	END DO
	2615	70 CONTINUE
[879]	2616
[1992]	2617	! -3/ height of triangle with area= sum_dth and base = dthmin
[879]	2618
[1992]	2619	hw = 2.*sum_dth/min(dthmin, -0.5)
	2620	hw = max(hwmin, hw)
[879]	2621
[1992]	2622	! -4/ now, get Ptop
[879]	2623
[1992]	2624	ktop = 0
	2625	z = 0.
[879]	2626
[1992]	2627	DO k = 1, klev
	2628	dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg), hw-z)
	2629	IF (dz<=0) GO TO 75
	2630	z = z + dz
	2631	ptop = ph(k) - rho(k)rgdz
	2632	ktop = k
	2633	END DO
	2634	75 CONTINUE
[879]	2635
[1992]	2636	! -5/Correct ktop and ptop
[879]	2637
[1992]	2638	ptop_new = ptop
[879]	2639
[1992]	2640	DO k = ktop, 2, -1
	2641	IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
	2642	ptop_new = ((dth(k)+delta_t_min)p(k-1)-(dth(k-1)+delta_t_min)p(k))/ &
	2643	(dth(k)-dth(k-1))
	2644	GO TO 275
	2645	END IF
	2646	END DO
	2647	275 CONTINUE
[879]	2648
[1992]	2649	ptop = ptop_new
[879]	2650
[1992]	2651	DO k = klev, 1, -1
	2652	IF (ph(k+1)<ptop) ktop = k
	2653	END DO
[879]	2654
[1992]	2655	! -6/ Set deltatw & deltaqw to 0 above kupper
[879]	2656
[1992]	2657	DO k = kupper, klev
	2658	deltatw(k) = 0.
	2659	deltaqw(k) = 0.
	2660	END DO
[879]	2661
[1992]	2662	! ------------------------------------------------------------------
	2663	END DO ! end sub-timestep loop
	2664	! -----------------------------------------------------------------
[879]	2665
[1992]	2666	! Get back to tendencies per second
[879]	2667
[1992]	2668	DO k = 1, kupper - 1
	2669	dtls(k) = dtls(k)/dtime
	2670	dqls(k) = dqls(k)/dtime
	2671	d_deltatw2(k) = d_deltatw2(k)/dtime
	2672	d_deltaqw2(k) = d_deltaqw2(k)/dtime
	2673	d_deltat_gw(k) = d_deltat_gw(k)/dtime
	2674	END DO
[879]	2675
[1992]	2676	! 2.1 - Undisturbed area and Wake integrals
	2677	! ---------------------------------------------------------
[879]	2678
[1992]	2679	z = 0.
	2680	sum_thu = 0.
	2681	sum_tu = 0.
	2682	sum_qu = 0.
	2683	sum_thvu = 0.
	2684	sum_dth = 0.
	2685	sum_dq = 0.
	2686	sum_rho = 0.
	2687	sum_dtdwn = 0.
	2688	sum_dqdwn = 0.
[879]	2689
[1992]	2690	av_thu = 0.
	2691	av_tu = 0.
	2692	av_qu = 0.
	2693	av_thvu = 0.
	2694	av_dth = 0.
	2695	av_dq = 0.
	2696	av_rho = 0.
	2697	av_dtdwn = 0.
	2698	av_dqdwn = 0.
[879]	2699
[1992]	2700	! Potential temperatures and humidity
	2701	! ----------------------------------------------------------
[879]	2702
[1992]	2703	DO k = 1, klev
	2704	rho(k) = p(k)/(rd*te(k))
	2705	IF (k==1) THEN
	2706	rhoh(k) = ph(k)/(rd*te(k))
	2707	zhh(k) = 0
	2708	ELSE
	2709	rhoh(k) = ph(k)2./(rd(te(k)+te(k-1)))
	2710	zhh(k) = (ph(k)-ph(k-1))/(-rhoh(k)*rg) + zhh(k-1)
	2711	END IF
	2712	the(k) = te(k)/ppi(k)
	2713	thu(k) = (te(k)-deltatw(k)*sigmaw)/ppi(k)
	2714	tu(k) = te(k) - deltatw(k)*sigmaw
	2715	qu(k) = qe(k) - deltaqw(k)*sigmaw
	2716	rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
	2717	dth(k) = deltatw(k)/ppi(k)
[879]	2718
[1992]	2719	END DO
[879]	2720
[1992]	2721	! Integrals (and wake top level number)
	2722	! -----------------------------------------------------------
[879]	2723
[1992]	2724	! Initialize sum_thvu to 1st level virt. pot. temp.
[879]	2725
[1992]	2726	z = 1.
	2727	dz = 1.
	2728	sum_thvu = thu(1)(1.+epsqu(1))*dz
	2729	sum_dth = 0.
[879]	2730
[1992]	2731	DO k = 1, klev
	2732	dz = -(max(ph(k+1),ptop)-ph(k))/(rho(k)*rg)
[879]	2733
[1992]	2734	IF (dz<=0) GO TO 51
	2735	z = z + dz
	2736	sum_thu = sum_thu + thu(k)*dz
	2737	sum_tu = sum_tu + tu(k)*dz
	2738	sum_qu = sum_qu + qu(k)*dz
	2739	sum_thvu = sum_thvu + thu(k)(1.+epsqu(k))*dz
	2740	sum_dth = sum_dth + dth(k)*dz
	2741	sum_dq = sum_dq + deltaqw(k)*dz
	2742	sum_rho = sum_rho + rhow(k)*dz
	2743	sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
	2744	sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
	2745	END DO
	2746	51 CONTINUE
[879]	2747
[1992]	2748	hw0 = z
[879]	2749
[1992]	2750	! 2.1 - WAPE and mean forcing computation
	2751	! -------------------------------------------------------------
[879]	2752
[1992]	2753	! Means
[879]	2754
[1992]	2755	av_thu = sum_thu/hw0
	2756	av_tu = sum_tu/hw0
	2757	av_qu = sum_qu/hw0
	2758	av_thvu = sum_thvu/hw0
	2759	av_dth = sum_dth/hw0
	2760	av_dq = sum_dq/hw0
	2761	av_rho = sum_rho/hw0
	2762	av_dtdwn = sum_dtdwn/hw0
	2763	av_dqdwn = sum_dqdwn/hw0
[879]	2764
[1992]	2765	wape2 = -rghw0(av_dth+eps(av_thuav_dq+av_dthav_qu+av_dthav_dq))/ &
	2766	av_thvu
[879]	2767
	2768
[1992]	2769	! 2.2 Prognostic variable update
	2770	! ------------------------------------------------------------
[879]	2771
[1992]	2772	! Filter out bad wakes
[879]	2773
[1992]	2774	IF (wape2<0.) THEN
	2775	IF (prt_level>=10) PRINT *, 'wape2<0'
	2776	wape2 = 0.
	2777	hw = hwmin
	2778	sigmaw = max(sigmad, sigd_con)
	2779	fip = 0.
	2780	DO k = 1, klev
	2781	deltatw(k) = 0.
	2782	deltaqw(k) = 0.
	2783	dth(k) = 0.
	2784	END DO
	2785	ELSE
	2786	IF (prt_level>=10) PRINT *, 'wape2>0'
	2787	cstar2 = starksqrt(2.wape2)
[879]	2788
[1992]	2789	END IF
[879]	2790
[1992]	2791	ktopw = ktop
[879]	2792
[1992]	2793	IF (ktopw>0) THEN
[879]	2794
[1992]	2795	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
	2796	! cc heff = 600.
	2797	! Utilisation de la hauteur hw
	2798	! c heff = 0.7*hw
	2799	heff = hw
[879]	2800
[1992]	2801	fip = 0.5rho(ktopw)cstar2*3heff2sqrt(sigmawwdens3.14)
	2802	fip = alpk*fip
	2803	! jyg2
	2804	ELSE
	2805	fip = 0.
	2806	END IF
[879]	2807
	2808
[1992]	2809	! Limitation de sigmaw
[879]	2810
[1992]	2811	! sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
	2812	! alors il disparait en se mélangeant à la partie undisturbed
[879]	2813
[1992]	2814	! correction NICOLAS . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
	2815	IF ((sigmaw>0.9) .OR. ((wape>=wape2) .AND. (wape2<= &
	2816	1.0)) .OR. (ktopw<=2)) THEN
	2817	! IM cf NR/JYG 251108 . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
	2818	! IF (sigmaw.GT.0.9) THEN
	2819	DO k = 1, klev
	2820	dtls(k) = 0.
	2821	dqls(k) = 0.
	2822	deltatw(k) = 0.
	2823	deltaqw(k) = 0.
	2824	END DO
	2825	wape = 0.
	2826	hw = hwmin
	2827	sigmaw = sigmad
	2828	fip = 0.
	2829	END IF
[879]	2830
[1992]	2831	RETURN
	2832	END SUBROUTINE wake_scal
[879]	2833
[1992]	2834
	2835

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

Download in other formats:

Original Format