Context Navigation

wake.F @ 1117

Last change on this file since 1117 was 1059, checked in by lmdzadmin, 16 years ago
Ajout condition supplementaire disparition poche froide (si elle n'atteint pas au - 2eme niveau) NR/JYG/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 67.0 KB

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

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

Download in other formats:

Original Format