Context Navigation

wake.F @ 1366

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

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

Download in other formats:

Original Format