Context Navigation

wake.F90 @ 2435

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

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

Download in other formats:

Original Format