Context Navigation

wake.F @ 2394

Last change on this file since 2394 was 1336, checked in by idelkadi, 15 years ago
Remise de la valeur par default de alp_offset a 0. Ajout de la possiblite de lecture des parametres de ini_wake.F dans le fichier ini_wake_param.data Ajout de la possiblite de lecture des parametres de wake.F dans le fichier wake_param.data Correction dans la partie convection (nouvelle physique)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 82.7 KB

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

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

Download in other formats:

Original Format