Context Navigation

wake.F90 @ 2093

Last change on this file since 2093 was 2078, checked in by idelkadi, 10 years ago
Rajout des directives OPENMP pour la lecture des cles et parametres de convection et wakes dans les fichier conv_param.data, ep_param.data et wake_param.data
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: 78.2 KB

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

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

Download in other formats:

Original Format