Context Navigation

wake.F @ 1155

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

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

Download in other formats:

Original Format