Context Navigation

wake.F90 @ 2686

Last change on this file since 2686 was 2671, checked in by jyg, 8 years ago
Bug fix in cv3_routine.F90 (cv3_unsat could create a precipitating downdraught when convection was off) and in wake.F90 and calwake.F90 (some array dimensions were incompatible and some initializations were missing)
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: 59.5 KB

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

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

Download in other formats:

Original Format