Context Navigation

wake.F90 @ 2157

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

Rev	Line
[1992]	1
[1403]	2	! $Id: wake.F90 1999 2014-03-20 09:57:19Z acaubel $
[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