Context Navigation

wake.F90 @ 2580

Last change on this file since 2580 was 2495, checked in by lguez, 9 years ago

Bug fix. Virtual temperature is:

[1 + (1 / eps - 1) q] T

not:

(1 + eps q) T

(no worry: the numerical difference is small).

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 53.0 KB

Rev	Line
[1992]	1
[1403]	2	! $Id: wake.F90 2495 2016-04-14 14:20:44Z idelkadi $
[879]	3
[2308]	4	SUBROUTINE wake(p, ph, pi, dtime, sigd_con, &
	5	te0, qe0, omgb, &
	6	dtdwn, dqdwn, amdwn, amup, dta, dqa, &
	7	wdtpbl, wdqpbl, udtpbl, udqpbl, &
	8	deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, &
	9	dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, &
	10	dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
	11	d_deltat_gw, d_deltatw2, d_deltaqw2)
[1146]	12
[974]	13
[1992]	14	! **************************************************************
	15	! *
	16	! WAKE *
	17	! retour a un Pupper fixe *
	18	! *
	19	! written by : GRANDPEIX Jean-Yves 09/03/2000 *
	20	! modified by : ROEHRIG Romain 01/29/2007 *
	21	! **************************************************************
[974]	22
[2474]	23	USE ioipsl_getin_p_mod, ONLY : getin_p
[1992]	24	USE dimphy
[2078]	25	use mod_phys_lmdz_para
[2311]	26	USE print_control_mod, ONLY: prt_level
[1992]	27	IMPLICIT NONE
	28	! ============================================================================
[974]	29
	30
[1992]	31	! But : Decrire le comportement des poches froides apparaissant dans les
	32	! grands systemes convectifs, et fournir l'energie disponible pour
	33	! le declenchement de nouvelles colonnes convectives.
[974]	34
[1992]	35	! Variables d'etat : deltatw : ecart de temperature wake-undisturbed
	36	! area
	37	! deltaqw : ecart d'humidite wake-undisturbed area
	38	! sigmaw : fraction d'aire occupee par la poche.
[974]	39
[1992]	40	! Variable de sortie :
[974]	41
[1992]	42	! wape : WAke Potential Energy
	43	! fip : Front Incident Power (W/m2) - ALP
	44	! gfl : Gust Front Length per unit area (m-1)
	45	! dtls : large scale temperature tendency due to wake
	46	! dqls : large scale humidity tendency due to wake
	47	! hw : hauteur de la poche
	48	! dp_omgb : vertical gradient of large scale omega
	49	! wdens : densite de poches
	50	! omgbdth: flux of Delta_Theta transported by LS omega
	51	! dtKE : differential heating (wake - unpertubed)
	52	! dqKE : differential moistening (wake - unpertubed)
	53	! omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
	54	! dp_deltomg : vertical gradient of omg (s-1)
	55	! spread : spreading term in dt_wake and dq_wake
	56	! deltatw : updated temperature difference (T_w-T_u).
	57	! deltaqw : updated humidity difference (q_w-q_u).
	58	! sigmaw : updated wake fractional area.
	59	! d_deltat_gw : delta T tendency due to GW
[974]	60
[1992]	61	! Variables d'entree :
[974]	62
[1992]	63	! aire : aire de la maille
	64	! te0 : temperature dans l'environnement (K)
	65	! qe0 : humidite dans l'environnement (kg/kg)
	66	! omgb : vitesse verticale moyenne sur la maille (Pa/s)
	67	! dtdwn: source de chaleur due aux descentes (K/s)
	68	! dqdwn: source d'humidite due aux descentes (kg/kg/s)
	69	! dta : source de chaleur due courants satures et detrain (K/s)
	70	! dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
	71	! amdwn: flux de masse total des descentes, par unite de
	72	! surface de la maille (kg/m2/s)
	73	! amup : flux de masse total des ascendances, par unite de
	74	! surface de la maille (kg/m2/s)
	75	! p : pressions aux milieux des couches (Pa)
	76	! ph : pressions aux interfaces (Pa)
	77	! pi : (p/p_0)**kapa (adim)
	78	! dtime: increment temporel (s)
[974]	79
[1992]	80	! Variables internes :
[974]	81
[1992]	82	! rhow : masse volumique de la poche froide
	83	! rho : environment density at P levels
	84	! rhoh : environment density at Ph levels
	85	! te : environment temperature \| may change within
	86	! qe : environment humidity \| sub-time-stepping
	87	! the : environment potential temperature
	88	! thu : potential temperature in undisturbed area
	89	! tu : temperature in undisturbed area
	90	! qu : humidity in undisturbed area
	91	! dp_omgb: vertical gradient og LS omega
	92	! omgbw : wake average vertical omega
	93	! dp_omgbw: vertical gradient of omgbw
	94	! omgbdq : flux of Delta_q transported by LS omega
	95	! dth : potential temperature diff. wake-undist.
	96	! th1 : first pot. temp. for vertical advection (=thu)
	97	! th2 : second pot. temp. for vertical advection (=thw)
	98	! q1 : first humidity for vertical advection
	99	! q2 : second humidity for vertical advection
	100	! d_deltatw : terme de redistribution pour deltatw
	101	! d_deltaqw : terme de redistribution pour deltaqw
	102	! deltatw0 : deltatw initial
	103	! deltaqw0 : deltaqw initial
	104	! hw0 : hw initial
	105	! sigmaw0: sigmaw initial
	106	! amflux : horizontal mass flux through wake boundary
	107	! wdens_ref: initial number of wakes per unit area (3D) or per
	108	! unit length (2D), at the beginning of each time step
	109	! Tgw : 1 sur la période de onde de gravité
	110	! Cgw : vitesse de propagation de onde de gravité
	111	! LL : distance entre 2 poches
[974]	112
[1992]	113	! -------------------------------------------------------------------------
	114	! Déclaration de variables
	115	! -------------------------------------------------------------------------
[1146]	116
[1992]	117	include "YOMCST.h"
	118	include "cvthermo.h"
[974]	119
[1992]	120	! Arguments en entree
	121	! --------------------
[974]	122
[2308]	123	REAL, DIMENSION (klon, klev), INTENT(IN) :: p, pi
	124	REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph, omgb
	125	REAL, INTENT(IN) :: dtime
	126	REAL, DIMENSION (klon, klev), INTENT(IN) :: te0, qe0
	127	REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn
	128	REAL, DIMENSION (klon, klev), INTENT(IN) :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED
	129	REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup
	130	REAL, DIMENSION (klon, klev), INTENT(IN) :: dta, dqa
	131	REAL, DIMENSION (klon), INTENT(IN) :: sigd_con
[974]	132
[2308]	133	!
	134	! Input/Output
	135	REAL, DIMENSION (klon, klev), INTENT(INOUT) :: deltatw, deltaqw
	136	REAL, DIMENSION (klon), INTENT(INOUT) :: sigmaw
	137
[1992]	138	! Sorties
	139	! --------
[974]	140
[2308]	141	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth
	142	REAL, DIMENSION (klon, klev), INTENT(OUT) :: tu, qu
	143	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls
	144	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke
	145	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtpbl, dqpbl
	146	REAL, DIMENSION (klon, klev), INTENT(OUT) :: spread
	147	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltatw2, d_deltaqw2
	148	REAL, DIMENSION (klon, klev+1), INTENT(OUT) :: omgbdth, omg
	149	REAL, DIMENSION (klon, klev), INTENT(OUT) :: dp_omgb, dp_deltomg
	150	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_gw
	151	REAL, DIMENSION (klon), INTENT(OUT) :: hw, wape, fip, gfl, cstar
	152	REAL, DIMENSION (klon), INTENT(OUT) :: wdens
	153	INTEGER, DIMENSION (klon), INTENT(OUT) :: ktopw
[974]	154
[1992]	155	! Variables internes
	156	! -------------------
[974]	157
[1992]	158	! Variables à fixer
	159	REAL alon
[2078]	160	LOGICAL, SAVE :: first = .TRUE.
	161	!$OMP THREADPRIVATE(first)
	162	REAL, SAVE :: stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol
	163	!$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
[2467]	164
[1992]	165	REAL delta_t_min
	166	INTEGER nsub
	167	REAL dtimesub
	168	REAL sigmad, hwmin, wapecut
	169	REAL :: sigmaw_max
	170	REAL :: dens_rate
	171	REAL wdens0
	172	! IM 080208
	173	LOGICAL, DIMENSION (klon) :: gwake
[974]	174
[1992]	175	! Variables de sauvegarde
	176	REAL, DIMENSION (klon, klev) :: deltatw0
	177	REAL, DIMENSION (klon, klev) :: deltaqw0
	178	REAL, DIMENSION (klon, klev) :: te, qe
	179	REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
[974]	180
[1992]	181	! Variables pour les GW
	182	REAL, DIMENSION (klon) :: ll
	183	REAL, DIMENSION (klon, klev) :: n2
	184	REAL, DIMENSION (klon, klev) :: cgw
	185	REAL, DIMENSION (klon, klev) :: tgw
[1403]	186
[1992]	187	! Variables liées au calcul de hw
	188	REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
	189	REAL, DIMENSION (klon) :: sum_dth
	190	REAL, DIMENSION (klon) :: dthmin
	191	REAL, DIMENSION (klon) :: z, dz, hw0
	192	INTEGER, DIMENSION (klon) :: ktop, kupper
[1403]	193
[1992]	194	! Sub-timestep tendencies and related variables
	195	REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
	196	REAL d_te(klon, klev), d_qe(klon, klev)
	197	REAL d_sigmaw(klon), alpha(klon)
	198	REAL q0_min(klon), q1_min(klon)
	199	LOGICAL wk_adv(klon), ok_qx_qw(klon)
	200	REAL epsilon
	201	DATA epsilon/1.E-15/
[974]	202
[1992]	203	! Autres variables internes
	204	INTEGER isubstep, k, i
[974]	205
[1992]	206	REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
	207	REAL, DIMENSION (klon) :: sum_dq, sum_rho
	208	REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
	209	REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
	210	REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
	211	REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
[974]	212
[1992]	213	REAL, DIMENSION (klon, klev) :: rho, rhow
	214	REAL, DIMENSION (klon, klev+1) :: rhoh
	215	REAL, DIMENSION (klon, klev) :: rhow_moyen
	216	REAL, DIMENSION (klon, klev) :: zh
	217	REAL, DIMENSION (klon, klev+1) :: zhh
	218	REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
[974]	219
[1992]	220	REAL, DIMENSION (klon, klev) :: the, thu
[974]	221
[1992]	222	! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
[974]	223
[1992]	224	REAL, DIMENSION (klon, klev+1) :: omgbw
	225	REAL, DIMENSION (klon) :: pupper
	226	REAL, DIMENSION (klon) :: omgtop
	227	REAL, DIMENSION (klon, klev) :: dp_omgbw
	228	REAL, DIMENSION (klon) :: ztop, dztop
	229	REAL, DIMENSION (klon, klev) :: alpha_up
[974]	230
[1992]	231	REAL, DIMENSION (klon) :: rre1, rre2
	232	REAL :: rrd1, rrd2
	233	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
	234	REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
	235	REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
	236	REAL, DIMENSION (klon, klev) :: omgbdq
[974]	237
[1992]	238	REAL, DIMENSION (klon) :: ff, gg
	239	REAL, DIMENSION (klon) :: wape2, cstar2, heff
[974]	240
[1992]	241	REAL, DIMENSION (klon, klev) :: crep
[974]	242
[1992]	243	REAL, DIMENSION (klon, klev) :: ppi
[974]	244
[1992]	245	! cc nrlmd
	246	REAL, DIMENSION (klon) :: death_rate, nat_rate
	247	REAL, DIMENSION (klon, klev) :: entr
	248	REAL, DIMENSION (klon, klev) :: detr
[974]	249
[1992]	250	! -------------------------------------------------------------------------
	251	! Initialisations
	252	! -------------------------------------------------------------------------
[974]	253
[1992]	254	! print*, 'wake initialisations'
[974]	255
[1992]	256	! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
	257	! -------------------------------------------------------------------------
[974]	258
[1992]	259	DATA wapecut, sigmad, hwmin/5., .02, 10./
	260	! cc nrlmd
	261	DATA sigmaw_max/0.4/
	262	DATA dens_rate/0.1/
	263	! cc
	264	! Longueur de maille (en m)
	265	! -------------------------------------------------------------------------
[974]	266
[1992]	267	! ALON = 3.e5
	268	alon = 1.E6
[974]	269
	270
[1992]	271	! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
[974]	272
[1992]	273	! coefgw : Coefficient pour les ondes de gravité
	274	! stark : Coefficient k dans Cstar=ksqrt(2WAPE)
	275	! wdens : Densité de poche froide par maille
	276	! -------------------------------------------------------------------------
[974]	277
[1992]	278	! cc nrlmd coefgw=10
	279	! coefgw=1
	280	! wdens0 = 1.0/(alon**2)
	281	! cc nrlmd wdens = 1.0/(alon**2)
	282	! cc nrlmd stark = 0.50
	283	! CRtest
	284	! cc nrlmd alpk=0.1
	285	! alpk = 1.0
	286	! alpk = 0.5
	287	! alpk = 0.05
[1146]	288
[2078]	289	if (first) then
[1992]	290	crep_upper = 0.9
	291	crep_sol = 1.0
[974]	292
[1992]	293	! cc nrlmd Lecture du fichier wake_param.data
[2474]	294	stark=0.33
	295	CALL getin_p('stark',stark)
	296	alpk=0.25
	297	CALL getin_p('alpk',alpk)
	298	wdens_ref=8.E-12
	299	CALL getin_p('wdens_ref',wdens_ref)
	300	coefgw=4.
	301	CALL getin_p('coefgw',coefgw)
	302
	303	WRITE(,) 'stark=', stark
	304	WRITE(,) 'alpk=', alpk
	305	WRITE(,) 'wdens_ref=', wdens_ref
	306	WRITE(,) 'coefgw=', coefgw
	307
[2078]	308	first=.false.
	309	endif
	310
[1992]	311	! Initialisation de toutes des densites a wdens_ref.
	312	! Les densites peuvent evoluer si les poches debordent
	313	! (voir au tout debut de la boucle sur les substeps)
	314	wdens = wdens_ref
[974]	315
[1992]	316	! print*,'stark',stark
	317	! print*,'alpk',alpk
	318	! print*,'wdens',wdens
	319	! print*,'coefgw',coefgw
	320	! cc
	321	! Minimum value for \|T_wake - T_undist\|. Used for wake top definition
	322	! -------------------------------------------------------------------------
[974]	323
[1992]	324	delta_t_min = 0.2
[974]	325
[1992]	326	! 1. - Save initial values and initialize tendencies
	327	! --------------------------------------------------
[974]	328
[1992]	329	DO k = 1, klev
	330	DO i = 1, klon
	331	ppi(i, k) = pi(i, k)
	332	deltatw0(i, k) = deltatw(i, k)
	333	deltaqw0(i, k) = deltaqw(i, k)
	334	te(i, k) = te0(i, k)
	335	qe(i, k) = qe0(i, k)
	336	dtls(i, k) = 0.
	337	dqls(i, k) = 0.
	338	d_deltat_gw(i, k) = 0.
	339	d_te(i, k) = 0.
	340	d_qe(i, k) = 0.
	341	d_deltatw(i, k) = 0.
	342	d_deltaqw(i, k) = 0.
	343	! IM 060508 beg
	344	d_deltatw2(i, k) = 0.
	345	d_deltaqw2(i, k) = 0.
	346	! IM 060508 end
	347	END DO
	348	END DO
	349	! sigmaw1=sigmaw
	350	! IF (sigd_con.GT.sigmaw1) THEN
	351	! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
	352	! ENDIF
	353	DO i = 1, klon
	354	! c sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
	355	sigmaw(i) = amax1(sigmaw(i), sigmad)
	356	sigmaw(i) = amin1(sigmaw(i), 0.99)
	357	sigmaw0(i) = sigmaw(i)
	358	wape(i) = 0.
	359	wape2(i) = 0.
	360	d_sigmaw(i) = 0.
	361	ktopw(i) = 0
	362	END DO
[974]	363
	364
[1992]	365	! 2. - Prognostic part
	366	! --------------------
[974]	367
	368
[1992]	369	! 2.1 - Undisturbed area and Wake integrals
	370	! ---------------------------------------------------------
[974]	371
[1992]	372	DO i = 1, klon
	373	z(i) = 0.
	374	ktop(i) = 0
	375	kupper(i) = 0
	376	sum_thu(i) = 0.
	377	sum_tu(i) = 0.
	378	sum_qu(i) = 0.
	379	sum_thvu(i) = 0.
	380	sum_dth(i) = 0.
	381	sum_dq(i) = 0.
	382	sum_rho(i) = 0.
	383	sum_dtdwn(i) = 0.
	384	sum_dqdwn(i) = 0.
[974]	385
[1992]	386	av_thu(i) = 0.
	387	av_tu(i) = 0.
	388	av_qu(i) = 0.
	389	av_thvu(i) = 0.
	390	av_dth(i) = 0.
	391	av_dq(i) = 0.
	392	av_rho(i) = 0.
	393	av_dtdwn(i) = 0.
	394	av_dqdwn(i) = 0.
	395	END DO
[974]	396
[1992]	397	! Distance between wakes
	398	DO i = 1, klon
	399	ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
	400	END DO
	401	! Potential temperatures and humidity
	402	! ----------------------------------------------------------
	403	DO k = 1, klev
	404	DO i = 1, klon
	405	! write(,)'wake 1',i,k,rd,te(i,k)
	406	rho(i, k) = p(i, k)/(rd*te(i,k))
	407	! write(,)'wake 2',rho(i,k)
	408	IF (k==1) THEN
	409	! write(,)'wake 3',i,k,rd,te(i,k)
	410	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
	411	! write(,)'wake 4',i,k,rd,te(i,k)
	412	zhh(i, k) = 0
	413	ELSE
	414	! write(,)'wake 5',rd,(te(i,k)+te(i,k-1))
	415	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
	416	! write(,)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
	417	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
	418	END IF
	419	! write(,)'wake 7',ppi(i,k)
	420	the(i, k) = te(i, k)/ppi(i, k)
	421	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
	422	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
	423	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
	424	! write(,)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
	425	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
	426	dth(i, k) = deltatw(i, k)/ppi(i, k)
	427	END DO
	428	END DO
[1403]	429
[1992]	430	DO k = 1, klev - 1
	431	DO i = 1, klon
	432	IF (k==1) THEN
	433	n2(i, k) = 0
	434	ELSE
	435	n2(i, k) = amax1(0., -rg*2/the(i,k)rho(i,k)*(the(i,k+1)-the(i, &
	436	k-1))/(p(i,k+1)-p(i,k-1)))
	437	END IF
	438	zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
[1403]	439
[1992]	440	cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
	441	tgw(i, k) = coefgw*cgw(i, k)/ll(i)
	442	END DO
	443	END DO
[974]	444
[1992]	445	DO i = 1, klon
	446	n2(i, klev) = 0
	447	zh(i, klev) = 0
	448	cgw(i, klev) = 0
	449	tgw(i, klev) = 0
	450	END DO
[974]	451
[1992]	452	! Calcul de la masse volumique moyenne de la colonne (bdlmd)
	453	! -----------------------------------------------------------------
[974]	454
[1992]	455	DO k = 1, klev
	456	DO i = 1, klon
	457	epaisseur1(i, k) = 0.
	458	epaisseur2(i, k) = 0.
	459	END DO
	460	END DO
[974]	461
[1992]	462	DO i = 1, klon
	463	epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
	464	epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
	465	rhow_moyen(i, 1) = rhow(i, 1)
	466	END DO
[974]	467
[1992]	468	DO k = 2, klev
	469	DO i = 1, klon
	470	epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
	471	epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
	472	rhow_moyen(i, k) = (rhow_moyen(i,k-1)epaisseur2(i,k-1)+rhow(i,k) &
	473	epaisseur1(i,k))/epaisseur2(i, k)
	474	END DO
	475	END DO
[974]	476
	477
[1992]	478	! Choose an integration bound well above wake top
	479	! -----------------------------------------------------------------
[974]	480
[1992]	481	! Pupper = 50000. ! melting level
	482	! Pupper = 60000.
	483	! Pupper = 80000. ! essais pour case_e
	484	DO i = 1, klon
	485	pupper(i) = 0.6*ph(i, 1)
	486	pupper(i) = max(pupper(i), 45000.)
	487	! cc Pupper(i) = 60000.
	488	END DO
[1146]	489
[1403]	490
[1992]	491	! Determine Wake top pressure (Ptop) from buoyancy integral
	492	! --------------------------------------------------------
[1403]	493
[1992]	494	! -1/ Pressure of the level where dth becomes less than delta_t_min.
	495
	496	DO i = 1, klon
	497	ptop_provis(i) = ph(i, 1)
	498	END DO
	499	DO k = 2, klev
	500	DO i = 1, klon
	501
	502	! IM v3JYG; ptop_provis(i).LT. ph(i,1)
	503
	504	IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
	505	ptop_provis(i)==ph(i,1)) THEN
	506	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	507	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	508	END IF
	509	END DO
	510	END DO
	511
	512	! -2/ dth integral
	513
	514	DO i = 1, klon
	515	sum_dth(i) = 0.
	516	dthmin(i) = -delta_t_min
	517	z(i) = 0.
	518	END DO
	519
	520	DO k = 1, klev
	521	DO i = 1, klon
	522	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
	523	IF (dz(i)>0) THEN
	524	z(i) = z(i) + dz(i)
	525	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	526	dthmin(i) = amin1(dthmin(i), dth(i,k))
	527	END IF
	528	END DO
	529	END DO
	530
	531	! -3/ height of triangle with area= sum_dth and base = dthmin
	532
	533	DO i = 1, klon
	534	hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
	535	hw0(i) = amax1(hwmin, hw0(i))
	536	END DO
	537
	538	! -4/ now, get Ptop
	539
	540	DO i = 1, klon
	541	z(i) = 0.
	542	ptop(i) = ph(i, 1)
	543	END DO
	544
	545	DO k = 1, klev
	546	DO i = 1, klon
	547	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
	548	IF (dz(i)>0) THEN
	549	z(i) = z(i) + dz(i)
	550	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
	551	END IF
	552	END DO
	553	END DO
	554
	555
	556	! -5/ Determination de ktop et kupper
	557
	558	DO k = klev, 1, -1
	559	DO i = 1, klon
	560	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	561	IF (ph(i,k+1)<pupper(i)) kupper(i) = k
	562	END DO
	563	END DO
	564
	565	! On evite kupper = 1 et kupper = klev
	566	DO i = 1, klon
	567	kupper(i) = max(kupper(i), 2)
	568	kupper(i) = min(kupper(i), klev-1)
	569	END DO
	570
	571
	572	! -6/ Correct ktop and ptop
	573
	574	DO i = 1, klon
	575	ptop_new(i) = ptop(i)
	576	END DO
	577	DO k = klev, 2, -1
	578	DO i = 1, klon
	579	IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
	580	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
	581	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	582	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	583	END IF
	584	END DO
	585	END DO
	586
	587	DO i = 1, klon
	588	ptop(i) = ptop_new(i)
	589	END DO
	590
	591	DO k = klev, 1, -1
	592	DO i = 1, klon
	593	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	594	END DO
	595	END DO
	596
	597	! -5/ Set deltatw & deltaqw to 0 above kupper
	598
	599	DO k = 1, klev
	600	DO i = 1, klon
	601	IF (k>=kupper(i)) THEN
	602	deltatw(i, k) = 0.
	603	deltaqw(i, k) = 0.
	604	END IF
	605	END DO
	606	END DO
	607
	608
	609	! Vertical gradient of LS omega
	610
	611	DO k = 1, klev
	612	DO i = 1, klon
	613	IF (k<=kupper(i)) THEN
	614	dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
	615	END IF
	616	END DO
	617	END DO
	618
	619	! Integrals (and wake top level number)
	620	! --------------------------------------
	621
	622	! Initialize sum_thvu to 1st level virt. pot. temp.
	623
	624	DO i = 1, klon
	625	z(i) = 1.
	626	dz(i) = 1.
[2495]	627	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
[1992]	628	sum_dth(i) = 0.
	629	END DO
	630
	631	DO k = 1, klev
	632	DO i = 1, klon
	633	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	634	IF (dz(i)>0) THEN
	635	z(i) = z(i) + dz(i)
	636	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	637	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	638	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
[2495]	639	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
[1992]	640	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	641	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	642	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	643	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	644	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	645	END IF
	646	END DO
	647	END DO
	648
	649	DO i = 1, klon
	650	hw0(i) = z(i)
	651	END DO
	652
	653
	654	! 2.1 - WAPE and mean forcing computation
	655	! ---------------------------------------
	656
	657	! ---------------------------------------
	658
	659	! Means
	660
	661	DO i = 1, klon
	662	av_thu(i) = sum_thu(i)/hw0(i)
	663	av_tu(i) = sum_tu(i)/hw0(i)
	664	av_qu(i) = sum_qu(i)/hw0(i)
	665	av_thvu(i) = sum_thvu(i)/hw0(i)
	666	! av_thve = sum_thve/hw0
	667	av_dth(i) = sum_dth(i)/hw0(i)
	668	av_dq(i) = sum_dq(i)/hw0(i)
	669	av_rho(i) = sum_rho(i)/hw0(i)
	670	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	671	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	672
[2495]	673	wape(i) = -rghw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i)+av_dth(i)*av_qu(i &
[1992]	674	)+av_dth(i)*av_dq(i)))/av_thvu(i)
	675	END DO
	676
	677	! 2.2 Prognostic variable update
	678	! ------------------------------
	679
	680	! Filter out bad wakes
	681
	682	DO k = 1, klev
	683	DO i = 1, klon
	684	IF (wape(i)<0.) THEN
	685	deltatw(i, k) = 0.
	686	deltaqw(i, k) = 0.
	687	dth(i, k) = 0.
	688	END IF
	689	END DO
	690	END DO
	691
	692	DO i = 1, klon
	693	IF (wape(i)<0.) THEN
	694	wape(i) = 0.
	695	cstar(i) = 0.
	696	hw(i) = hwmin
	697	sigmaw(i) = amax1(sigmad, sigd_con(i))
	698	fip(i) = 0.
	699	gwake(i) = .FALSE.
	700	ELSE
	701	cstar(i) = starksqrt(2.wape(i))
	702	gwake(i) = .TRUE.
	703	END IF
	704	END DO
	705
	706
	707	! Check qx and qw positivity
	708	! --------------------------
	709	DO i = 1, klon
	710	q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, &
	711	1)+(1.-sigmaw(i))*deltaqw(i,1)))
	712	END DO
	713	DO k = 2, klev
	714	DO i = 1, klon
	715	q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, &
	716	k)+(1.-sigmaw(i))*deltaqw(i,k)))
	717	IF (q1_min(i)<=q0_min(i)) THEN
	718	q0_min(i) = q1_min(i)
	719	END IF
	720	END DO
	721	END DO
	722
	723	DO i = 1, klon
	724	ok_qx_qw(i) = q0_min(i) >= 0.
	725	alpha(i) = 1.
	726	END DO
	727
	728	! C -----------------------------------------------------------------
	729	! Sub-time-stepping
	730	! -----------------
	731
	732	nsub = 10
	733	dtimesub = dtime/nsub
	734
	735	! ------------------------------------------------------------
	736	DO isubstep = 1, nsub
	737	! ------------------------------------------------------------
	738
	739	! wk_adv is the logical flag enabling wake evolution in the time advance
	740	! loop
	741	DO i = 1, klon
	742	wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
	743	END DO
	744
	745	! cc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement
	746	! négatif de ktop à kupper --------
	747	! cc On calcule pour cela une densité wdens0 pour laquelle on
	748	! aurait un entrainement nul ---
	749	DO i = 1, klon
	750	! c print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
	751	! c $ isubstep,wk_adv(i),cstar(i),wape(i)
	752	IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
	753	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
	754	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
	755	wdens0 = (sigmaw(i)/(4.3.14))((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( &
	756	ph(i,1)-pupper(i))cstar(i)))*(2)
	757	IF (wdens(i)<=wdens0*1.1) THEN
	758	wdens(i) = wdens0
	759	END IF
	760	! c print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
	761	! c $ ,ph(i,1)-pupper(i)',
	762	! c $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
	763	! c $ ,ph(i,1)-pupper(i)
	764	END IF
	765	END DO
	766
	767	! cc nrlmd
	768
	769	DO i = 1, klon
	770	IF (wk_adv(i)) THEN
[1403]	771	gfl(i) = 2.sqrt(3.14wdens(i)*sigmaw(i))
[1992]	772	sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
	773	END IF
	774	END DO
	775	DO i = 1, klon
	776	IF (wk_adv(i)) THEN
	777	! cc nrlmd Introduction du taux de mortalité des poches et
	778	! test sur sigmaw_max=0.4
	779	! cc d_sigmaw(i) = gfl(i)Cstar(i)dtimesub
	780	IF (sigmaw(i)>=sigmaw_max) THEN
	781	death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
	782	ELSE
	783	death_rate(i) = 0.
	784	END IF
	785	d_sigmaw(i) = gfl(i)cstar(i)dtimesub - death_rate(i)sigmaw(i) &
	786	dtimesub
	787	! $ - nat_rate(i)sigmaw(i)dtimesub
	788	! c print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
	789	! c $ death_rate(i),ktop(i),kupper(i)',
	790	! c $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
	791	! c $ death_rate(i),ktop(i),kupper(i)
[1403]	792
[1992]	793	! sigmaw(i) =sigmaw(i) + gfl(i)Cstar(i)dtimesub
	794	! sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
	795	! wdens = wdens0/(10.*sigmaw)
	796	! sigmaw =max(sigmaw,sigd_con)
	797	! sigmaw =max(sigmaw,sigmad)
	798	END IF
	799	END DO
	800
	801
	802	! calcul de la difference de vitesse verticale poche - zone non perturbee
	803	! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
	804	! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on
	805	! definit
	806	! IM 060208 au niveau k=1..?
	807	DO k = 1, klev
	808	DO i = 1, klon
	809	IF (wk_adv(i)) THEN !!! nrlmd
	810	dp_deltomg(i, k) = 0.
	811	END IF
	812	END DO
	813	END DO
	814	DO k = 1, klev + 1
	815	DO i = 1, klon
	816	IF (wk_adv(i)) THEN !!! nrlmd
	817	omg(i, k) = 0.
	818	END IF
	819	END DO
	820	END DO
	821
	822	DO i = 1, klon
	823	IF (wk_adv(i)) THEN
	824	z(i) = 0.
	825	omg(i, 1) = 0.
	826	dp_deltomg(i, 1) = -(gfl(i)cstar(i))/(sigmaw(i)(1-sigmaw(i)))
	827	END IF
	828	END DO
	829
	830	DO k = 2, klev
	831	DO i = 1, klon
	832	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
[974]	833	dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
[1992]	834	z(i) = z(i) + dz(i)
	835	dp_deltomg(i, k) = dp_deltomg(i, 1)
	836	omg(i, k) = dp_deltomg(i, 1)*z(i)
	837	END IF
	838	END DO
	839	END DO
	840
	841	DO i = 1, klon
	842	IF (wk_adv(i)) THEN
	843	dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
	844	ztop(i) = z(i) + dztop(i)
	845	omgtop(i) = dp_deltomg(i, 1)*ztop(i)
	846	END IF
	847	END DO
	848
	849	! -----------------
	850	! From m/s to Pa/s
	851	! -----------------
	852
	853	DO i = 1, klon
	854	IF (wk_adv(i)) THEN
	855	omgtop(i) = -rho(i, ktop(i))rgomgtop(i)
	856	dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
	857	END IF
	858	END DO
	859
	860	DO k = 1, klev
	861	DO i = 1, klon
	862	IF (wk_adv(i) .AND. k<=ktop(i)) THEN
	863	omg(i, k) = -rho(i, k)rgomg(i, k)
	864	dp_deltomg(i, k) = dp_deltomg(i, 1)
	865	END IF
	866	END DO
	867	END DO
	868
	869	! raccordement lineaire de omg de ptop a pupper
	870
	871	DO i = 1, klon
	872	IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
	873	omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
	874	rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
	875	dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
	876	(ptop(i)-pupper(i))
	877	END IF
	878	END DO
	879
	880	! c DO i=1,klon
	881	! c print*,'Pente entre 0 et kupper (référence)'
	882	! c $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
	883	! c print*,'Pente entre ktop et kupper'
	884	! c $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
	885	! c ENDDO
	886	! c
	887	DO k = 1, klev
	888	DO i = 1, klon
	889	IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
	890	dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
	891	omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
	892	END IF
	893	END DO
	894	END DO
	895	! cc nrlmd
	896	! c DO i=1,klon
	897	! c print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
	898	! c END DO
	899	! cc
	900
	901
	902	! -- Compute wake average vertical velocity omgbw
	903
	904
	905	DO k = 1, klev + 1
	906	DO i = 1, klon
[1146]	907	IF (wk_adv(i)) THEN
[1992]	908	omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
	909	END IF
	910	END DO
	911	END DO
	912	! -- and its vertical gradient dp_omgbw
	913
	914	DO k = 1, klev
	915	DO i = 1, klon
[1146]	916	IF (wk_adv(i)) THEN
[1992]	917	dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
	918	END IF
	919	END DO
	920	END DO
[974]	921
[1992]	922	! -- Upstream coefficients for omgb velocity
	923	! -- (alpha_up(k) is the coefficient of the value at level k)
	924	! -- (1-alpha_up(k) is the coefficient of the value at level k-1)
	925	DO k = 1, klev
	926	DO i = 1, klon
	927	IF (wk_adv(i)) THEN
	928	alpha_up(i, k) = 0.
	929	IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
	930	END IF
	931	END DO
	932	END DO
[974]	933
[1992]	934	! Matrix expressing [The,deltatw] from [Th1,Th2]
[974]	935
[1992]	936	DO i = 1, klon
	937	IF (wk_adv(i)) THEN
	938	rre1(i) = 1. - sigmaw(i)
	939	rre2(i) = sigmaw(i)
	940	END IF
	941	END DO
	942	rrd1 = -1.
	943	rrd2 = 1.
[974]	944
[1992]	945	! -- Get [Th1,Th2], dth and [q1,q2]
[974]	946
[1992]	947	DO k = 1, klev
	948	DO i = 1, klon
	949	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
	950	dth(i, k) = deltatw(i, k)/ppi(i, k)
	951	th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
	952	th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
	953	q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
	954	q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
	955	END IF
	956	END DO
	957	END DO
[974]	958
[1992]	959	DO i = 1, klon
	960	IF (wk_adv(i)) THEN !!! nrlmd
	961	d_th1(i, 1) = 0.
	962	d_th2(i, 1) = 0.
	963	d_dth(i, 1) = 0.
	964	d_q1(i, 1) = 0.
	965	d_q2(i, 1) = 0.
	966	d_dq(i, 1) = 0.
	967	END IF
	968	END DO
[974]	969
[1992]	970	DO k = 2, klev
	971	DO i = 1, klon
	972	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
	973	d_th1(i, k) = th1(i, k-1) - th1(i, k)
	974	d_th2(i, k) = th2(i, k-1) - th2(i, k)
	975	d_dth(i, k) = dth(i, k-1) - dth(i, k)
	976	d_q1(i, k) = q1(i, k-1) - q1(i, k)
	977	d_q2(i, k) = q2(i, k-1) - q2(i, k)
	978	d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
	979	END IF
	980	END DO
	981	END DO
[1146]	982
[1992]	983	DO i = 1, klon
	984	IF (wk_adv(i)) THEN
	985	omgbdth(i, 1) = 0.
	986	omgbdq(i, 1) = 0.
	987	END IF
	988	END DO
[1277]	989
[1992]	990	DO k = 2, klev
	991	DO i = 1, klon
	992	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN ! loop on interfaces
	993	omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
	994	omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
	995	END IF
	996	END DO
	997	END DO
[1403]	998
[1992]	999	! -----------------------------------------------------------------
	1000	DO k = 1, klev
	1001	DO i = 1, klon
	1002	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
	1003	! -----------------------------------------------------------------
[974]	1004
[1992]	1005	! Compute redistribution (advective) term
[1403]	1006
[1992]	1007	d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
	1008	(rrd1omg(i,k)sigmaw(i)d_th1(i,k)-rrd2omg(i,k+1)*(1.-sigmaw( &
	1009	i))d_th2(i,k+1)-(1.-alpha_up(i,k))omgbdth(i,k)-alpha_up(i,k+1)* &
	1010	omgbdth(i,k+1))*ppi(i, k)
	1011	! print*,'d_deltatw=',d_deltatw(i,k)
[1403]	1012
[1992]	1013	d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
	1014	(rrd1omg(i,k)sigmaw(i)d_q1(i,k)-rrd2omg(i,k+1)*(1.-sigmaw( &
	1015	i))d_q2(i,k+1)-(1.-alpha_up(i,k))omgbdq(i,k)-alpha_up(i,k+1)* &
	1016	omgbdq(i,k+1))
	1017	! print*,'d_deltaqw=',d_deltaqw(i,k)
[974]	1018
[1992]	1019	! and increment large scale tendencies
[974]	1020
	1021
	1022
	1023
[1992]	1024	! C
	1025	! -----------------------------------------------------------------
	1026	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i, &
	1027	k)-rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
	1028	k+1)) & ! cc nrlmd $
	1029	! -sigmaw(i)(1.-sigmaw(i))dth(i,k)*dp_deltomg(i,k)
	1030	-sigmaw(i)(1.-sigmaw(i))dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
	1031	k)-ph(i,k+1)) & ! cc
	1032	)*ppi(i, k)
[974]	1033
[1992]	1034	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i, &
	1035	k)-rre2(i)omg(i,k+1)(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
	1036	k+1)) & ! cc nrlmd $
	1037	! -sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*dp_deltomg(i,k)
	1038	-sigmaw(i)(1.-sigmaw(i))deltaqw(i,k)*(omg(i,k)-omg(i, &
	1039	k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
	1040	)
	1041	! cc nrlmd
	1042	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
	1043	d_te(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_th1(i, &
	1044	k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k)
[1403]	1045
[1992]	1046	d_qe(i, k) = dtimesub((rre1(i)omg(i,k)sigmaw(i)d_q1(i, &
	1047	k)/(ph(i,k)-ph(i,k+1))))
[1403]	1048
[1992]	1049	END IF
	1050	! cc
	1051	END DO
	1052	END DO
	1053	! ------------------------------------------------------------------
[974]	1054
[1992]	1055	! Increment state variables
[974]	1056
[1992]	1057	DO k = 1, klev
	1058	DO i = 1, klon
	1059	! cc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
	1060	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1061	! cc
[974]	1062
[1146]	1063
[974]	1064
[1992]	1065	! Coefficient de répartition
[974]	1066
[1992]	1067	crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
	1068	(ph(i,kupper(i))-ph(i,1))
	1069	crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i &
	1070	,kupper(i)))
[974]	1071
	1072
[1992]	1073	! Reintroduce compensating subsidence term.
[1146]	1074
[1992]	1075	! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
	1076	! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
	1077	! . /(1-sigmaw)
	1078	! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
	1079	! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
	1080	! . /(1-sigmaw)
[974]	1081
[1992]	1082	! dtKE(k)=(dtdwn(k)Crep(k)+(1-Crep(k))dta(k))/sigmaw
	1083	! dtKE(k)=dtKE(k)-(dtdwn(k)(1-Crep(k))+dta(k)Crep(k))
	1084	! . /(1-sigmaw)
	1085	! dqKE(k)=(dqdwn(k)Crep(k)+(1-Crep(k))dqa(k))/sigmaw
	1086	! dqKE(k)=dqKE(k)-(dqdwn(k)(1-Crep(k))+dqa(k)Crep(k))
	1087	! . /(1-sigmaw)
[974]	1088
[1992]	1089	dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
	1090	dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
	1091	! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
[974]	1092
[2155]	1093	!jyg<
	1094	!!
	1095	!!---------------------------------------------------------------
	1096	!! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
	1097	!! is accounted for by the PBL and the Thermals schemes. It is now set to zero
	1098	!! within the Wake scheme.
	1099	!!---------------------------------------------------------------
	1100	dtPBL(i,k) = 0.
	1101	dqPBL(i,k) = 0.
	1102	!
	1103	!! dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
	1104	!! dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
	1105	!
	1106	!! dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
	1107	!! dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
[1992]	1108	! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
[2155]	1109	!>jyg
	1110	!
[1146]	1111
[1992]	1112	! cc nrlmd Prise en compte du taux de mortalité
	1113	! cc Définitions de entr, detr
	1114	detr(i, k) = 0.
[1146]	1115
[1992]	1116	entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
	1117	sigmaw(i)(1.-sigmaw(i))dp_deltomg(i, k)
[1146]	1118
[1992]	1119	spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
	1120	! cc spread(i,k) =
	1121	! (1.-sigmaw(i))dp_deltomg(i,k)+gfl(i)Cstar(i)/
	1122	! cc $ sigmaw(i)
[1146]	1123
	1124
[1992]	1125	! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
	1126	! Jingmei
[1146]	1127
[1992]	1128	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
	1129	! & Tgw(i,k),deltatw(i,k)
	1130	d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)deltatw(i, k) &
	1131	dtimesub
	1132	! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
	1133	ff(i) = d_deltatw(i, k)/dtimesub
[1403]	1134
[1992]	1135	! Sans GW
[1403]	1136
[1992]	1137	! deltatw(k)=deltatw(k)+dtimesub(ff+dtKE(k)-spread(k)deltatw(k))
[974]	1138
[1992]	1139	! GW formule 1
	1140
	1141	! deltatw(k) = deltatw(k)+dtimesub*
	1142	! $ (ff+dtKE(k) - spread(k)deltatw(k)-Tgw(k)deltatw(k))
	1143
	1144	! GW formule 2
	1145
	1146	IF (dtimesub*tgw(i,k)<1.E-10) THEN
	1147	d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
	1148	! $
	1149	! -spread(i,k)*deltatw(i,k)
	1150	-entr(i,k)deltatw(i,k)/sigmaw(i)-(death_rate(i)sigmaw( &
	1151	i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
	1152	-tgw(i,k)*deltatw(i,k))
	1153	ELSE
	1154	d_deltatw(i, k) = 1/tgw(i, k)(1-exp(-dtimesubtgw(i, &
	1155	k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc $
	1156	! -spread(i,k)*deltatw(i,k)
	1157	-entr(i,k)deltatw(i,k)/sigmaw(i)-(death_rate(i)sigmaw( &
	1158	i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
	1159	-tgw(i,k)*deltatw(i,k))
	1160	END IF
	1161
	1162	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1163
	1164	gg(i) = d_deltaqw(i, k)/dtimesub
	1165
	1166	d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc $
	1167	! -spread(i,k)*deltaqw(i,k))
	1168	-entr(i,k)deltaqw(i,k)/sigmaw(i)-(death_rate(i)sigmaw(i)+detr( &
	1169	i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
	1170	! cc
	1171
	1172	! cc nrlmd
	1173	! cc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
	1174	! cc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
	1175	! cc
	1176	END IF
	1177	END DO
	1178	END DO
	1179
	1180
	1181	! Scale tendencies so that water vapour remains positive in w and x.
	1182
	1183	CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
	1184	d_deltaqw, sigmaw, d_sigmaw, alpha)
	1185
	1186	! cc nrlmd
	1187	! c print*,'alpha'
	1188	! c do i=1,klon
	1189	! c print*,alpha(i)
	1190	! c end do
	1191	! cc
	1192	DO k = 1, klev
	1193	DO i = 1, klon
	1194	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1195	d_te(i, k) = alpha(i)*d_te(i, k)
	1196	d_qe(i, k) = alpha(i)*d_qe(i, k)
	1197	d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
	1198	d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
	1199	d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
	1200	END IF
	1201	END DO
	1202	END DO
	1203	DO i = 1, klon
	1204	IF (wk_adv(i)) THEN
	1205	d_sigmaw(i) = alpha(i)*d_sigmaw(i)
	1206	END IF
	1207	END DO
	1208
	1209	! Update large scale variables and wake variables
	1210	! IM 060208 manque DO i + remplace DO k=1,kupper(i)
	1211	! IM 060208 DO k = 1,kupper(i)
	1212	DO k = 1, klev
	1213	DO i = 1, klon
	1214	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1215	dtls(i, k) = dtls(i, k) + d_te(i, k)
	1216	dqls(i, k) = dqls(i, k) + d_qe(i, k)
	1217	! cc nrlmd
	1218	d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
	1219	d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
	1220	! cc
	1221	END IF
	1222	END DO
	1223	END DO
	1224	DO k = 1, klev
	1225	DO i = 1, klon
	1226	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
	1227	te(i, k) = te0(i, k) + dtls(i, k)
	1228	qe(i, k) = qe0(i, k) + dqls(i, k)
	1229	the(i, k) = te(i, k)/ppi(i, k)
	1230	deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
	1231	deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
	1232	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1233	! c print,'k,qx,qw',k,qe(i,k)-sigmaw(i)deltaqw(i,k)
	1234	! c $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
	1235	END IF
	1236	END DO
	1237	END DO
	1238	DO i = 1, klon
	1239	IF (wk_adv(i)) THEN
	1240	sigmaw(i) = sigmaw(i) + d_sigmaw(i)
	1241	END IF
	1242	END DO
	1243
	1244
	1245	! Determine Ptop from buoyancy integral
	1246	! ---------------------------------------
	1247
	1248	! - 1/ Pressure of the level where dth changes sign.
	1249
	1250	DO i = 1, klon
	1251	IF (wk_adv(i)) THEN
	1252	ptop_provis(i) = ph(i, 1)
	1253	END IF
	1254	END DO
	1255
	1256	DO k = 2, klev
	1257	DO i = 1, klon
	1258	IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
	1259	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
	1260	ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	1261	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	1262	END IF
	1263	END DO
	1264	END DO
	1265
	1266	! - 2/ dth integral
	1267
	1268	DO i = 1, klon
	1269	IF (wk_adv(i)) THEN !!! nrlmd
	1270	sum_dth(i) = 0.
	1271	dthmin(i) = -delta_t_min
[974]	1272	z(i) = 0.
[1992]	1273	END IF
	1274	END DO
	1275
	1276	DO k = 1, klev
	1277	DO i = 1, klon
	1278	IF (wk_adv(i)) THEN
	1279	dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
	1280	IF (dz(i)>0) THEN
	1281	z(i) = z(i) + dz(i)
	1282	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1283	dthmin(i) = amin1(dthmin(i), dth(i,k))
	1284	END IF
	1285	END IF
	1286	END DO
	1287	END DO
	1288
	1289	! - 3/ height of triangle with area= sum_dth and base = dthmin
	1290
	1291	DO i = 1, klon
	1292	IF (wk_adv(i)) THEN
	1293	hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
	1294	hw(i) = amax1(hwmin, hw(i))
	1295	END IF
	1296	END DO
	1297
	1298	! - 4/ now, get Ptop
	1299
	1300	DO i = 1, klon
	1301	IF (wk_adv(i)) THEN !!! nrlmd
	1302	ktop(i) = 0
	1303	z(i) = 0.
	1304	END IF
	1305	END DO
	1306
	1307	DO k = 1, klev
	1308	DO i = 1, klon
	1309	IF (wk_adv(i)) THEN
	1310	dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
	1311	IF (dz(i)>0) THEN
	1312	z(i) = z(i) + dz(i)
	1313	ptop(i) = ph(i, k) - rho(i, k)rgdz(i)
	1314	ktop(i) = k
	1315	END IF
	1316	END IF
	1317	END DO
	1318	END DO
	1319
	1320	! 4.5/Correct ktop and ptop
	1321
	1322	DO i = 1, klon
	1323	IF (wk_adv(i)) THEN
	1324	ptop_new(i) = ptop(i)
	1325	END IF
	1326	END DO
	1327
	1328	DO k = klev, 2, -1
	1329	DO i = 1, klon
	1330	! IM v3JYG; IF (k .GE. ktop(i)
	1331	IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
	1332	dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
	1333	ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
	1334	k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
	1335	END IF
	1336	END DO
	1337	END DO
	1338
	1339
	1340	DO i = 1, klon
	1341	IF (wk_adv(i)) THEN
	1342	ptop(i) = ptop_new(i)
	1343	END IF
	1344	END DO
	1345
	1346	DO k = klev, 1, -1
	1347	DO i = 1, klon
	1348	IF (wk_adv(i)) THEN !!! nrlmd
	1349	IF (ph(i,k+1)<ptop(i)) ktop(i) = k
	1350	END IF
	1351	END DO
	1352	END DO
	1353
	1354	! 5/ Set deltatw & deltaqw to 0 above kupper
	1355
	1356	DO k = 1, klev
	1357	DO i = 1, klon
	1358	IF (wk_adv(i) .AND. k>=kupper(i)) THEN
	1359	deltatw(i, k) = 0.
	1360	deltaqw(i, k) = 0.
	1361	END IF
	1362	END DO
	1363	END DO
	1364
	1365
	1366	! -------------Cstar computation---------------------------------
	1367	DO i = 1, klon
	1368	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1369	sum_thu(i) = 0.
	1370	sum_tu(i) = 0.
	1371	sum_qu(i) = 0.
	1372	sum_thvu(i) = 0.
	1373	sum_dth(i) = 0.
	1374	sum_dq(i) = 0.
	1375	sum_rho(i) = 0.
	1376	sum_dtdwn(i) = 0.
	1377	sum_dqdwn(i) = 0.
	1378
	1379	av_thu(i) = 0.
[1992]	1380	av_tu(i) = 0.
	1381	av_qu(i) = 0.
[974]	1382	av_thvu(i) = 0.
	1383	av_dth(i) = 0.
	1384	av_dq(i) = 0.
[1992]	1385	av_rho(i) = 0.
	1386	av_dtdwn(i) = 0.
[974]	1387	av_dqdwn(i) = 0.
[1992]	1388	END IF
	1389	END DO
[974]	1390
[1992]	1391	! Integrals (and wake top level number)
	1392	! --------------------------------------
[974]	1393
[1992]	1394	! Initialize sum_thvu to 1st level virt. pot. temp.
[974]	1395
[1992]	1396	DO i = 1, klon
	1397	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1398	z(i) = 1.
	1399	dz(i) = 1.
[2495]	1400	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
[974]	1401	sum_dth(i) = 0.
[1992]	1402	END IF
	1403	END DO
[974]	1404
[1992]	1405	DO k = 1, klev
	1406	DO i = 1, klon
	1407	IF (wk_adv(i)) THEN !!! nrlmd
	1408	dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	1409	IF (dz(i)>0) THEN
	1410	z(i) = z(i) + dz(i)
	1411	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	1412	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	1413	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
[2495]	1414	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
[1992]	1415	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1416	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	1417	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	1418	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	1419	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	1420	END IF
	1421	END IF
	1422	END DO
	1423	END DO
	1424
	1425	DO i = 1, klon
	1426	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1427	hw0(i) = z(i)
[1992]	1428	END IF
	1429	END DO
[974]	1430
	1431
[1992]	1432	! - WAPE and mean forcing computation
	1433	! ---------------------------------------
	1434
	1435	! ---------------------------------------
	1436
	1437	! Means
	1438
	1439	DO i = 1, klon
	1440	IF (wk_adv(i)) THEN !!! nrlmd
[974]	1441	av_thu(i) = sum_thu(i)/hw0(i)
	1442	av_tu(i) = sum_tu(i)/hw0(i)
	1443	av_qu(i) = sum_qu(i)/hw0(i)
	1444	av_thvu(i) = sum_thvu(i)/hw0(i)
	1445	av_dth(i) = sum_dth(i)/hw0(i)
	1446	av_dq(i) = sum_dq(i)/hw0(i)
	1447	av_rho(i) = sum_rho(i)/hw0(i)
	1448	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	1449	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	1450
[2495]	1451	wape(i) = -rghw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i)+av_dth(i)* &
[1992]	1452	av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
	1453	END IF
	1454	END DO
[974]	1455
[1992]	1456	! Filter out bad wakes
[974]	1457
[1992]	1458	DO k = 1, klev
	1459	DO i = 1, klon
	1460	IF (wk_adv(i)) THEN !!! nrlmd
	1461	IF (wape(i)<0.) THEN
	1462	deltatw(i, k) = 0.
	1463	deltaqw(i, k) = 0.
	1464	dth(i, k) = 0.
	1465	END IF
	1466	END IF
	1467	END DO
	1468	END DO
[974]	1469
[1992]	1470	DO i = 1, klon
	1471	IF (wk_adv(i)) THEN !!! nrlmd
	1472	IF (wape(i)<0.) THEN
	1473	wape(i) = 0.
	1474	cstar(i) = 0.
	1475	hw(i) = hwmin
	1476	sigmaw(i) = max(sigmad, sigd_con(i))
	1477	fip(i) = 0.
	1478	gwake(i) = .FALSE.
	1479	ELSE
	1480	cstar(i) = starksqrt(2.wape(i))
	1481	gwake(i) = .TRUE.
	1482	END IF
	1483	END IF
	1484	END DO
	1485
	1486	END DO ! end sub-timestep loop
	1487
	1488	! -----------------------------------------------------------------
	1489	! Get back to tendencies per second
	1490
	1491	DO k = 1, klev
	1492	DO i = 1, klon
	1493
	1494	! cc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
	1495	IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
	1496	! cc
	1497	dtls(i, k) = dtls(i, k)/dtime
	1498	dqls(i, k) = dqls(i, k)/dtime
	1499	d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
	1500	d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
	1501	d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
	1502	! c print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
	1503	! c $ ,death_rate(i)*sigmaw(i)
	1504	END IF
	1505	END DO
	1506	END DO
	1507
	1508
	1509	! ----------------------------------------------------------
	1510	! Determine wake final state; recompute wape, cstar, ktop;
	1511	! filter out bad wakes.
	1512	! ----------------------------------------------------------
	1513
	1514	! 2.1 - Undisturbed area and Wake integrals
	1515	! ---------------------------------------------------------
	1516
	1517	DO i = 1, klon
	1518	! cc nrlmd if (wk_adv(i)) then !!! nrlmd
	1519	IF (ok_qx_qw(i)) THEN
	1520	! cc
	1521	z(i) = 0.
	1522	sum_thu(i) = 0.
	1523	sum_tu(i) = 0.
	1524	sum_qu(i) = 0.
	1525	sum_thvu(i) = 0.
	1526	sum_dth(i) = 0.
	1527	sum_dq(i) = 0.
	1528	sum_rho(i) = 0.
	1529	sum_dtdwn(i) = 0.
	1530	sum_dqdwn(i) = 0.
	1531
	1532	av_thu(i) = 0.
	1533	av_tu(i) = 0.
	1534	av_qu(i) = 0.
	1535	av_thvu(i) = 0.
	1536	av_dth(i) = 0.
	1537	av_dq(i) = 0.
	1538	av_rho(i) = 0.
	1539	av_dtdwn(i) = 0.
	1540	av_dqdwn(i) = 0.
	1541	END IF
	1542	END DO
	1543	! Potential temperatures and humidity
	1544	! ----------------------------------------------------------
	1545
	1546	DO k = 1, klev
	1547	DO i = 1, klon
	1548	! cc nrlmd IF ( wk_adv(i)) THEN
	1549	IF (ok_qx_qw(i)) THEN
	1550	! cc
	1551	rho(i, k) = p(i, k)/(rd*te(i,k))
	1552	IF (k==1) THEN
	1553	rhoh(i, k) = ph(i, k)/(rd*te(i,k))
	1554	zhh(i, k) = 0
	1555	ELSE
	1556	rhoh(i, k) = ph(i, k)2./(rd(te(i,k)+te(i,k-1)))
	1557	zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
	1558	END IF
	1559	the(i, k) = te(i, k)/ppi(i, k)
	1560	thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
	1561	tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
	1562	qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
	1563	rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
	1564	dth(i, k) = deltatw(i, k)/ppi(i, k)
	1565	END IF
	1566	END DO
	1567	END DO
	1568
	1569	! Integrals (and wake top level number)
	1570	! -----------------------------------------------------------
	1571
	1572	! Initialize sum_thvu to 1st level virt. pot. temp.
	1573
	1574	DO i = 1, klon
	1575	! cc nrlmd IF ( wk_adv(i)) THEN
	1576	IF (ok_qx_qw(i)) THEN
	1577	! cc
	1578	z(i) = 1.
	1579	dz(i) = 1.
[2495]	1580	sum_thvu(i) = thu(i, 1)(1.+epsim1qu(i,1))*dz(i)
[1992]	1581	sum_dth(i) = 0.
	1582	END IF
	1583	END DO
	1584
	1585	DO k = 1, klev
	1586	DO i = 1, klon
	1587	! cc nrlmd IF ( wk_adv(i)) THEN
	1588	IF (ok_qx_qw(i)) THEN
	1589	! cc
	1590	dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
	1591	IF (dz(i)>0) THEN
	1592	z(i) = z(i) + dz(i)
	1593	sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
	1594	sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
	1595	sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
[2495]	1596	sum_thvu(i) = sum_thvu(i) + thu(i, k)(1.+epsim1qu(i,k))*dz(i)
[1992]	1597	sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
	1598	sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
	1599	sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
	1600	sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
	1601	sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
	1602	END IF
	1603	END IF
	1604	END DO
	1605	END DO
	1606
	1607	DO i = 1, klon
	1608	! cc nrlmd IF ( wk_adv(i)) THEN
	1609	IF (ok_qx_qw(i)) THEN
	1610	! cc
	1611	hw0(i) = z(i)
	1612	END IF
	1613	END DO
	1614
	1615	! - WAPE and mean forcing computation
	1616	! -------------------------------------------------------------
	1617
	1618	! Means
	1619
	1620	DO i = 1, klon
	1621	! cc nrlmd IF ( wk_adv(i)) THEN
	1622	IF (ok_qx_qw(i)) THEN
	1623	! cc
	1624	av_thu(i) = sum_thu(i)/hw0(i)
	1625	av_tu(i) = sum_tu(i)/hw0(i)
	1626	av_qu(i) = sum_qu(i)/hw0(i)
	1627	av_thvu(i) = sum_thvu(i)/hw0(i)
	1628	av_dth(i) = sum_dth(i)/hw0(i)
	1629	av_dq(i) = sum_dq(i)/hw0(i)
	1630	av_rho(i) = sum_rho(i)/hw0(i)
	1631	av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
	1632	av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
	1633
[2495]	1634	wape2(i) = -rghw0(i)(av_dth(i)+epsim1(av_thu(i)av_dq(i)+av_dth(i)* &
[1992]	1635	av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
	1636	END IF
	1637	END DO
	1638
	1639	! Prognostic variable update
	1640	! ------------------------------------------------------------
	1641
	1642	! Filter out bad wakes
	1643
	1644	DO k = 1, klev
	1645	DO i = 1, klon
	1646	! cc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
	1647	IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
	1648	! cc
	1649	deltatw(i, k) = 0.
	1650	deltaqw(i, k) = 0.
	1651	dth(i, k) = 0.
	1652	END IF
	1653	END DO
	1654	END DO
	1655
	1656
	1657	DO i = 1, klon
	1658	! cc nrlmd IF ( wk_adv(i)) THEN
	1659	IF (ok_qx_qw(i)) THEN
	1660	! cc
	1661	IF (wape2(i)<0.) THEN
[974]	1662	wape2(i) = 0.
[1992]	1663	cstar2(i) = 0.
[974]	1664	hw(i) = hwmin
[1992]	1665	sigmaw(i) = amax1(sigmad, sigd_con(i))
[974]	1666	fip(i) = 0.
	1667	gwake(i) = .FALSE.
	1668	ELSE
[1992]	1669	IF (prt_level>=10) PRINT *, 'wape2>0'
	1670	cstar2(i) = starksqrt(2.wape2(i))
[974]	1671	gwake(i) = .TRUE.
[1992]	1672	END IF
	1673	END IF
	1674	END DO
[974]	1675
[1992]	1676	DO i = 1, klon
	1677	! cc nrlmd IF ( wk_adv(i)) THEN
	1678	IF (ok_qx_qw(i)) THEN
	1679	! cc
	1680	ktopw(i) = ktop(i)
	1681	END IF
	1682	END DO
[974]	1683
[1992]	1684	DO i = 1, klon
	1685	! cc nrlmd IF ( wk_adv(i)) THEN
	1686	IF (ok_qx_qw(i)) THEN
	1687	! cc
	1688	IF (ktopw(i)>0 .AND. gwake(i)) THEN
[1403]	1689
[1992]	1690	! jyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
	1691	! cc heff = 600.
	1692	! Utilisation de la hauteur hw
	1693	! c heff = 0.7*hw
	1694	heff(i) = hw(i)
[1403]	1695
[1992]	1696	fip(i) = 0.5rho(i, ktopw(i))cstar2(i)*3heff(i)2 &
	1697	sqrt(sigmaw(i)wdens(i)3.14)
	1698	fip(i) = alpk*fip(i)
	1699	! jyg2
	1700	ELSE
	1701	fip(i) = 0.
	1702	END IF
	1703	END IF
	1704	END DO
[1146]	1705
[1992]	1706	! Limitation de sigmaw
	1707
	1708	! cc nrlmd
	1709	! DO i=1,klon
	1710	! IF (OK_qx_qw(i)) THEN
	1711	! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
	1712	! ENDIF
	1713	! ENDDO
	1714	! cc
	1715	DO k = 1, klev
	1716	DO i = 1, klon
	1717
	1718	! cc nrlmd On maintient désormais constant sigmaw en régime
	1719	! permanent
	1720	! cc IF ((sigmaw(i).GT.sigmaw_max).or.
	1721	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
	1722	.NOT. ok_qx_qw(i)) THEN
	1723	! cc
	1724	dtls(i, k) = 0.
	1725	dqls(i, k) = 0.
	1726	deltatw(i, k) = 0.
	1727	deltaqw(i, k) = 0.
	1728	END IF
	1729	END DO
	1730	END DO
	1731
	1732	! cc nrlmd On maintient désormais constant sigmaw en régime permanent
	1733	DO i = 1, klon
	1734	IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
	1735	.NOT. ok_qx_qw(i)) THEN
	1736	wape(i) = 0.
	1737	cstar(i) = 0.
[2308]	1738	!!jyg Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
	1739	!! hw(i) = hwmin !jyg
	1740	!! sigmaw(i) = sigmad !jyg
	1741	hw(i) = 0. !jyg
	1742	sigmaw(i) = 0. !jyg
[1992]	1743	fip(i) = 0.
	1744	ELSE
	1745	wape(i) = wape2(i)
	1746	cstar(i) = cstar2(i)
	1747	END IF
	1748	! c print*,'wape wape2 ktopw OK_qx_qw =',
	1749	! c $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
	1750	END DO
	1751
	1752
	1753	RETURN
	1754	END SUBROUTINE wake
	1755
	1756	SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
	1757	d_deltaqw, sigmaw, d_sigmaw, alpha)
	1758	! ------------------------------------------------------
	1759	! Dtermination du coefficient alpha tel que les tendances
	1760	! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
	1761	! a une humidite positive dans la zone (x) et dans la zone (w).
	1762	! ------------------------------------------------------
[2197]	1763	IMPLICIT NONE
[1992]	1764
	1765	! Input
	1766	REAL qe(nlon, nl), d_qe(nlon, nl)
	1767	REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
	1768	REAL sigmaw(nlon), d_sigmaw(nlon)
	1769	LOGICAL wk_adv(nlon)
	1770	INTEGER nl, nlon
	1771	! Output
	1772	REAL alpha(nlon)
	1773	! Internal variables
	1774	REAL zeta(nlon, nl)
	1775	REAL alpha1(nlon)
	1776	REAL x, a, b, c, discrim
	1777	REAL epsilon
	1778	! DATA epsilon/1.e-15/
[2197]	1779	INTEGER i,k
[1992]	1780
	1781	DO k = 1, nl
	1782	DO i = 1, nlon
	1783	IF (wk_adv(i)) THEN
	1784	IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
	1785	zeta(i, k) = 0.
[1146]	1786	ELSE
[1992]	1787	zeta(i, k) = 1.
[1146]	1788	END IF
[1992]	1789	END IF
	1790	END DO
	1791	DO i = 1, nlon
	1792	IF (wk_adv(i)) THEN
	1793	x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
	1794	(zeta(i,k)-sigmaw(i))d_deltaqw(i, k) - d_sigmaw(i)(deltaqw(i,k)+ &
	1795	d_deltaqw(i,k))
	1796	a = -d_sigmaw(i)*d_deltaqw(i, k)
	1797	b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
	1798	deltaqw(i, k)*d_sigmaw(i)
	1799	c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
	1800	discrim = bb - 4.a*c
	1801	! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
	1802	IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
	1803	alpha1(i) = 1.
[1146]	1804	ELSE
[1992]	1805	IF (x>=0.) THEN
	1806	alpha1(i) = 1.
	1807	ELSE
	1808	IF (a>0.) THEN
	1809	alpha1(i) = 0.9min((2.c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
	1810	))/(2.*a))
	1811	ELSE IF (a==0.) THEN
	1812	alpha1(i) = 0.9*(-c/b)
	1813	ELSE
	1814	! print*,'a,b,c discrim',a,b,c discrim
	1815	alpha1(i) = 0.9max((2.c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
	1816	))/(2.*a))
	1817	END IF
	1818	END IF
	1819	END IF
	1820	alpha(i) = min(alpha(i), alpha1(i))
	1821	END IF
	1822	END DO
	1823	END DO
[1146]	1824
[1992]	1825	RETURN
	1826	END SUBROUTINE wake_vec_modulation
[974]	1827
[879]	1828

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

Download in other formats:

Original Format