Context Navigation

orografi.F90 @ 5403

Last change on this file since 5403 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: 43.7 KB

Rev	Line
[1992]	1
[1403]	2	! $Id: orografi.F90 2408 2015-12-14 10:43:09Z abarral $
[524]	3
[1992]	4	SUBROUTINE drag_noro(nlon, nlev, dtime, paprs, pplay, pmea, pstd, psig, pgam, &
	5	pthe, ppic, pval, kgwd, kdx, ktest, t, u, v, pulow, pvlow, pustr, pvstr, &
	6	d_t, d_u, d_v)
[524]	7
[1992]	8	USE dimphy
	9	IMPLICIT NONE
	10	! ======================================================================
	11	! Auteur(s): F.Lott (LMD/CNRS) date: 19950201
	12	! Objet: Frottement de la montagne Interface
	13	! ======================================================================
	14	! Arguments:
	15	! dtime---input-R- pas d'integration (s)
	16	! paprs---input-R-pression pour chaque inter-couche (en Pa)
	17	! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	18	! t-------input-R-temperature (K)
	19	! u-------input-R-vitesse horizontale (m/s)
	20	! v-------input-R-vitesse horizontale (m/s)
[524]	21
[1992]	22	! d_t-----output-R-increment de la temperature
	23	! d_u-----output-R-increment de la vitesse u
	24	! d_v-----output-R-increment de la vitesse v
	25	! ======================================================================
	26	include "YOMCST.h"
[524]	27
[1992]	28	! ARGUMENTS
[524]	29
[1992]	30	INTEGER nlon, nlev
	31	REAL dtime
	32	REAL paprs(klon, klev+1)
	33	REAL pplay(klon, klev)
	34	REAL pmea(nlon), pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon)
	35	REAL ppic(nlon), pval(nlon)
	36	REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon)
	37	REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev)
	38	REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev)
[524]	39
[1992]	40	INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
[524]	41
[1992]	42	! Variables locales:
[524]	43
[1992]	44	REAL zgeom(klon, klev)
	45	REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev)
	46	REAL pt(klon, klev), pu(klon, klev), pv(klon, klev)
	47	REAL papmf(klon, klev), papmh(klon, klev+1)
[524]	48
[1992]	49	! initialiser les variables de sortie (pour securite)
[524]	50
[1992]	51	DO i = 1, klon
	52	pulow(i) = 0.0
	53	pvlow(i) = 0.0
	54	pustr(i) = 0.0
	55	pvstr(i) = 0.0
	56	END DO
	57	DO k = 1, klev
	58	DO i = 1, klon
	59	d_t(i, k) = 0.0
	60	d_u(i, k) = 0.0
	61	d_v(i, k) = 0.0
	62	pdudt(i, k) = 0.0
	63	pdvdt(i, k) = 0.0
	64	pdtdt(i, k) = 0.0
	65	END DO
	66	END DO
[524]	67
[1992]	68	! preparer les variables d'entree (attention: l'ordre des niveaux
	69	! verticaux augmente du haut vers le bas)
[524]	70
[1992]	71	DO k = 1, klev
	72	DO i = 1, klon
	73	pt(i, k) = t(i, klev-k+1)
	74	pu(i, k) = u(i, klev-k+1)
	75	pv(i, k) = v(i, klev-k+1)
	76	papmf(i, k) = pplay(i, klev-k+1)
	77	END DO
	78	END DO
	79	DO k = 1, klev + 1
	80	DO i = 1, klon
	81	papmh(i, k) = paprs(i, klev-k+2)
	82	END DO
	83	END DO
	84	DO i = 1, klon
	85	zgeom(i, klev) = rdpt(i, klev)log(papmh(i,klev+1)/papmf(i,klev))
	86	END DO
	87	DO k = klev - 1, 1, -1
	88	DO i = 1, klon
	89	zgeom(i, k) = zgeom(i, k+1) + rd(pt(i,k)+pt(i,k+1))/2.0log(papmf(i,k+ &
	90	1)/papmf(i,k))
	91	END DO
	92	END DO
[524]	93
[1992]	94	! appeler la routine principale
[524]	95
[1992]	96	CALL orodrag(klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, zgeom, pt, &
	97	pu, pv, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, pvlow, pdudt, &
	98	pdvdt, pdtdt)
[524]	99
[1992]	100	DO k = 1, klev
	101	DO i = 1, klon
	102	d_u(i, klev+1-k) = dtime*pdudt(i, k)
	103	d_v(i, klev+1-k) = dtime*pdvdt(i, k)
	104	d_t(i, klev+1-k) = dtime*pdtdt(i, k)
	105	pustr(i) = pustr(i) & ! IM BUG .
	106	! +rgpdudt(i,k)(papmh(i,k+1)-papmh(i,k))
	107	+pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
	108	pvstr(i) = pvstr(i) & ! IM BUG .
	109	! +rgpdvdt(i,k)(papmh(i,k+1)-papmh(i,k))
	110	+pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
	111	END DO
	112	END DO
[524]	113
[1992]	114	RETURN
	115	END SUBROUTINE drag_noro
	116	SUBROUTINE orodrag(nlon, nlev, kgwd, kdx, ktest, ptsphy, paphm1, papm1, &
	117	pgeom1, ptm1, pum1, pvm1, pmea, pstd, psig, pgamma, ptheta, ppic, pval &
	118	! outputs
	119	, pulow, pvlow, pvom, pvol, pte)
[524]	120
[1992]	121	USE dimphy
	122	IMPLICIT NONE
[524]	123
	124
	125
[1992]	126	! **** gwdrag - does the gravity wave parametrization.
[524]	127
[1992]	128	! purpose.
	129	! --------
[524]	130
[1992]	131	! this routine computes the physical tendencies of the
	132	! prognostic variables u,v and t due to vertical transports by
	133	! subgridscale orographically excited gravity waves
[524]	134
[1992]	135	! ** interface.
	136	! ----------
	137	! called from callpar.
[524]	138
[1992]	139	! the routine takes its input from the long-term storage:
	140	! u,v,t and p at t-1.
[524]	141
[1992]	142	! explicit arguments :
	143	! --------------------
	144	! ==== inputs ===
	145	! ==== outputs ===
[524]	146
[1992]	147	! implicit arguments : none
	148	! --------------------
[524]	149
[1992]	150	! implicit logical (l)
[524]	151
[1992]	152	! method.
	153	! -------
[524]	154
[1992]	155	! externals.
	156	! ----------
	157	INTEGER ismin, ismax
	158	EXTERNAL ismin, ismax
[524]	159
[1992]	160	! reference.
	161	! ----------
[524]	162
[1992]	163	! author.
	164	! -------
	165	! m.miller + b.ritter e.c.m.w.f. 15/06/86.
[524]	166
[1992]	167	! f.lott + m. miller e.c.m.w.f. 22/11/94
	168	! -----------------------------------------------------------------------
[524]	169
	170
[1992]	171	include "YOMCST.h"
	172	include "YOEGWD.h"
	173	! -----------------------------------------------------------------------
[524]	174
[1992]	175	! * 0.1 arguments
	176	! ---------
[524]	177
	178
[1992]	179	! ym integer nlon, nlev, klevm1
	180	INTEGER nlon, nlev
	181	INTEGER kgwd, jl, ilevp1, jk, ji
	182	REAL zdelp, ztemp, zforc, ztend
	183	REAL rover, zb, zc, zconb, zabsv
	184	REAL zzd1, ratio, zbet, zust, zvst, zdis
	185	REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(klon), &
	186	pvlow(klon)
	187	REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), pmea(nlon), &
	188	pstd(nlon), psig(nlon), pgamma(nlon), ptheta(nlon), ppic(nlon), &
	189	pval(nlon), pgeom1(nlon, nlev), papm1(nlon, nlev), paphm1(nlon, nlev+1)
[524]	190
[1992]	191	INTEGER kdx(nlon), ktest(nlon)
	192	! -----------------------------------------------------------------------
[524]	193
[1992]	194	! * 0.2 local arrays
	195	! ------------
	196	INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), iknu(klon), &
	197	iknu2(klon), ikcrit(klon), ikhlim(klon)
[524]	198
[1992]	199	REAL ztau(klon, klev+1), ztauf(klon, klev+1), zstab(klon, klev+1), &
	200	zvph(klon, klev+1), zrho(klon, klev+1), zri(klon, klev+1), &
	201	zpsi(klon, klev+1), zzdep(klon, klev)
	202	REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), &
	203	znu(klon), zd1(klon), zd2(klon), zdmod(klon)
	204	REAL ztmst, ptsphy, zrtmst
[524]	205
[1992]	206	! ------------------------------------------------------------------
[524]	207
[1992]	208	! * 1. initialization
	209	! --------------
[524]	210
	211
[1992]	212	! ------------------------------------------------------------------
[524]	213
[1992]	214	! * 1.1 computational constants
	215	! -----------------------
[524]	216
	217
[1992]	218	! ztmst=twodt
	219	! if(nstep.eq.nstart) ztmst=0.5*twodt
	220	! ym klevm1=klev-1
	221	ztmst = ptsphy
	222	zrtmst = 1./ztmst
[524]	223
[1992]	224	! ------------------------------------------------------------------
[524]	225
[1992]	226	! * 1.3 check whether row contains point for printing
	227	! ---------------------------------------------
[524]	228
	229
[1992]	230	! ------------------------------------------------------------------
[524]	231
[1992]	232	! * 2. precompute basic state variables.
	233	! * ---------- ----- ----- ----------
	234	! * define low level wind, project winds in plane of
	235	! * low level wind, determine sector in which to take
	236	! * the variance and set indicator for critical levels.
[524]	237
	238
	239
	240
[1992]	241	CALL orosetup(nlon, ktest, ikcrit, ikcrith, icrit, ikenvh, iknu, iknu2, &
	242	paphm1, papm1, pum1, pvm1, ptm1, pgeom1, pstd, zrho, zri, zstab, ztau, &
	243	zvph, zpsi, zzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, znu, &
	244	zd1, zd2, zdmod)
[524]	245
[1992]	246	! ***********************************************************
[524]	247
	248
[1992]	249	! * 3. compute low level stresses using subcritical and
	250	! * supercritical forms.computes anisotropy coefficient
	251	! * as measure of orographic twodimensionality.
[524]	252
	253
[1992]	254	CALL gwstress(nlon, nlev, ktest, icrit, ikenvh, iknu, zrho, zstab, zvph, &
	255	pstd, psig, pmea, ppic, ztau, pgeom1, zdmod)
	256
	257	! * 4. compute stress profile.
	258	! * ------- ------ --------
	259
	260
	261
	262	CALL gwprofil(nlon, nlev, kgwd, kdx, ktest, ikcrith, icrit, paphm1, zrho, &
	263	zstab, zvph, zri, ztau, zdmod, psig, pstd)
	264
	265	! * 5. compute tendencies.
	266	! * -------------------
	267
	268
	269	! explicit solution at all levels for the gravity wave
	270	! implicit solution for the blocked levels
	271
	272	DO jl = kidia, kfdia
	273	zvidis(jl) = 0.0
	274	zdudt(jl) = 0.0
	275	zdvdt(jl) = 0.0
	276	zdtdt(jl) = 0.0
	277	END DO
	278
	279	ilevp1 = klev + 1
	280
	281
	282	DO jk = 1, klev
	283
	284
	285	! do 523 jl=1,kgwd
	286	! ji=kdx(jl)
	287	! Modif vectorisation 02/04/2004
	288	DO ji = kidia, kfdia
	289	IF (ktest(ji)==1) THEN
	290
	291	zdelp = paphm1(ji, jk+1) - paphm1(ji, jk)
	292	ztemp = -rg(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)zdelp)
	293	zdudt(ji) = (pulow(ji)zd1(ji)-pvlow(ji)zd2(ji))*ztemp/zdmod(ji)
	294	zdvdt(ji) = (pvlow(ji)zd1(ji)+pulow(ji)zd2(ji))*ztemp/zdmod(ji)
	295
	296	! controle des overshoots:
	297
	298	zforc = sqrt(zdudt(ji)2+zdvdt(ji)2) + 1.E-12
	299	ztend = sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/ztmst + 1.E-12
	300	rover = 0.25
	301	IF (zforc>=rover*ztend) THEN
	302	zdudt(ji) = roverztend/zforczdudt(ji)
	303	zdvdt(ji) = roverztend/zforczdvdt(ji)
	304	END IF
	305
	306	! fin du controle des overshoots
	307
	308	IF (jk>=ikenvh(ji)) THEN
	309	zb = 1.0 - 0.18pgamma(ji) - 0.04pgamma(ji)**2
	310	zc = 0.48pgamma(ji) + 0.3pgamma(ji)**2
	311	zconb = 2.ztmstgkwakepsig(ji)/(4.pstd(ji))
	312	zabsv = sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/2.
	313	zzd1 = zbcos(zpsi(ji,jk))2 + zcsin(zpsi(ji,jk))**2
	314	ratio = (cos(zpsi(ji,jk))*2+pgamma(ji)sin(zpsi(ji, &
	315	jk))*2)/(pgamma(ji)cos(zpsi(ji,jk))2+sin(zpsi(ji,jk))2)
	316	zbet = max(0., 2.-1./ratio)zconbzzdep(ji, jk)zzd1zabsv
	317
	318	! simplement oppose au vent
	319
	320	zdudt(ji) = -pum1(ji, jk)/ztmst
	321	zdvdt(ji) = -pvm1(ji, jk)/ztmst
	322
	323	! projection dans la direction de l'axe principal de l'orographie
	324	! mod zdudt(ji)=-(pum1(ji,jk)cos(ptheta(ji)rpi/180.)
	325	! mod * +pvm1(ji,jk)sin(ptheta(ji)rpi/180.))
	326	! mod * cos(ptheta(ji)rpi/180.)/ztmst
	327	! mod zdvdt(ji)=-(pum1(ji,jk)cos(ptheta(ji)rpi/180.)
	328	! mod * +pvm1(ji,jk)sin(ptheta(ji)rpi/180.))
	329	! mod * sin(ptheta(ji)rpi/180.)/ztmst
	330	zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
	331	zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
	332	END IF
	333	pvom(ji, jk) = zdudt(ji)
	334	pvol(ji, jk) = zdvdt(ji)
	335	zust = pum1(ji, jk) + ztmst*zdudt(ji)
	336	zvst = pvm1(ji, jk) + ztmst*zdvdt(ji)
	337	zdis = 0.5(pum1(ji,jk)2+pvm1(ji,jk)2-zust2-zvst*2)
	338	zdedt(ji) = zdis/ztmst
	339	zvidis(ji) = zvidis(ji) + zdis*zdelp
	340	zdtdt(ji) = zdedt(ji)/rcpd
	341	! pte(ji,jk)=zdtdt(ji)
	342
	343	! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
	344
	345	pte(ji, jk) = 0.0
	346
	347	END IF
	348	END DO
	349
	350	END DO
	351
	352
	353	RETURN
	354	END SUBROUTINE orodrag
	355	SUBROUTINE orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, kkenvh, kknu, kknu2, &
	356	paphm1, papm1, pum1, pvm1, ptm1, pgeom1, pstd, prho, pri, pstab, ptau, &
	357	pvph, ppsi, pzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, pnu, &
	358	pd1, pd2, pdmod)
	359
	360	! **** gwsetup
	361
	362	! purpose.
	363	! --------
	364
	365	! ** interface.
	366	! ----------
	367	! from orodrag
	368
	369	! explicit arguments :
	370	! --------------------
	371	! ==== inputs ===
	372	! ==== outputs ===
	373
	374	! implicit arguments : none
	375	! --------------------
	376
	377	! method.
	378	! -------
	379
	380
	381	! externals.
	382	! ----------
	383
	384
	385	! reference.
	386	! ----------
	387
	388	! see ecmwf research department documentation of the "i.f.s."
	389
	390	! author.
	391	! -------
	392
	393	! modifications.
	394	! --------------
	395	! f.lott for the new-gwdrag scheme november 1993
	396
	397	! -----------------------------------------------------------------------
	398	USE dimphy
	399	IMPLICIT NONE
	400
	401
	402	include "YOMCST.h"
	403	include "YOEGWD.h"
	404
	405	! -----------------------------------------------------------------------
	406
	407	! * 0.1 arguments
	408	! ---------
	409
	410	INTEGER nlon
	411	INTEGER jl, jk
	412	REAL zdelp
	413
	414	INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ktest(nlon), kkenvh(nlon)
	415
	416
	417	REAL paphm1(nlon, klev+1), papm1(nlon, klev), pum1(nlon, klev), &
	418	pvm1(nlon, klev), ptm1(nlon, klev), pgeom1(nlon, klev), &
	419	prho(nlon, klev+1), pri(nlon, klev+1), pstab(nlon, klev+1), &
	420	ptau(nlon, klev+1), pvph(nlon, klev+1), ppsi(nlon, klev+1), &
	421	pzdep(nlon, klev)
	422	REAL pulow(nlon), pvlow(nlon), ptheta(nlon), pgamma(nlon), pnu(nlon), &
	423	pd1(nlon), pd2(nlon), pdmod(nlon)
	424	REAL pstd(nlon), pmea(nlon), ppic(nlon), pval(nlon)
	425
	426	! -----------------------------------------------------------------------
	427
	428	! * 0.2 local arrays
	429	! ------------
	430
	431
	432	INTEGER ilevm1, ilevm2, ilevh
	433	REAL zcons1, zcons2, zcons3, zhgeo
	434	REAL zu, zphi, zvt1, zvt2, zst, zvar, zdwind, zwind
	435	REAL zstabm, zstabp, zrhom, zrhop, alpha
	436	REAL zggeenv, zggeom1, zgvar
	437	LOGICAL lo
	438	LOGICAL ll1(klon, klev+1)
	439	INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), &
	440	ncount(klon)
	441
	442	REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev)
	443	REAL znorm(klon), zb(klon), zc(klon), zulow(klon), zvlow(klon), znup(klon), &
	444	znum(klon)
	445
	446	! ------------------------------------------------------------------
	447
	448	! * 1. initialization
	449	! --------------
	450
	451	! print *,' entree gwsetup'
	452
	453	! ------------------------------------------------------------------
	454
	455	! * 1.1 computational constants
	456	! -----------------------
	457
	458
	459	ilevm1 = klev - 1
	460	ilevm2 = klev - 2
	461	ilevh = klev/3
	462
	463	zcons1 = 1./rd
	464	! old zcons2=g**2/cpd
	465	zcons2 = rg**2/rcpd
	466	! old zcons3=1.5*api
	467	zcons3 = 1.5*rpi
	468
	469	! ------------------------------------------------------------------
	470
	471	! * 2.
	472	! --------------
	473
	474
	475	! ------------------------------------------------------------------
	476
	477	! * 2.1 define low level wind, project winds in plane of
	478	! * low level wind, determine sector in which to take
	479	! * the variance and set indicator for critical levels.
	480
	481
	482
	483	DO jl = kidia, kfdia
	484	kknu(jl) = klev
	485	kknu2(jl) = klev
	486	kknub(jl) = klev
	487	kknul(jl) = klev
	488	pgamma(jl) = max(pgamma(jl), gtsec)
	489	ll1(jl, klev+1) = .FALSE.
	490	END DO
	491
	492	! Ajouter une initialisation (L. Li, le 23fev99):
	493
	494	DO jk = klev, ilevh, -1
	495	DO jl = kidia, kfdia
	496	ll1(jl, jk) = .FALSE.
	497	END DO
	498	END DO
	499
	500	! * define top of low level flow
	501	! ----------------------------
	502	DO jk = klev, ilevh, -1
	503	DO jl = kidia, kfdia
	504	lo = (paphm1(jl,jk)/paphm1(jl,klev+1)) >= gsigcr
	505	IF (lo) THEN
	506	kkcrit(jl) = jk
	507	END IF
	508	zhcrit(jl, jk) = ppic(jl)
	509	zhgeo = pgeom1(jl, jk)/rg
	510	ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
	511	IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
	512	kknu(jl) = jk
	513	END IF
	514	IF (.NOT. ll1(jl,ilevh)) kknu(jl) = ilevh
	515	END DO
	516	END DO
	517	DO jk = klev, ilevh, -1
	518	DO jl = kidia, kfdia
	519	zhcrit(jl, jk) = ppic(jl) - pval(jl)
	520	zhgeo = pgeom1(jl, jk)/rg
	521	ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
	522	IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
	523	kknu2(jl) = jk
	524	END IF
	525	IF (.NOT. ll1(jl,ilevh)) kknu2(jl) = ilevh
	526	END DO
	527	END DO
	528	DO jk = klev, ilevh, -1
	529	DO jl = kidia, kfdia
	530	zhcrit(jl, jk) = amax1(ppic(jl)-pmea(jl), pmea(jl)-pval(jl))
	531	zhgeo = pgeom1(jl, jk)/rg
	532	ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
	533	IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
	534	kknub(jl) = jk
	535	END IF
	536	IF (.NOT. ll1(jl,ilevh)) kknub(jl) = ilevh
	537	END DO
	538	END DO
	539
	540	DO jl = kidia, kfdia
	541	kknu(jl) = min(kknu(jl), nktopg)
	542	kknu2(jl) = min(kknu2(jl), nktopg)
	543	kknub(jl) = min(kknub(jl), nktopg)
	544	kknul(jl) = klev
	545	END DO
	546
	547	! c* initialize various arrays
	548
	549	DO jl = kidia, kfdia
	550	prho(jl, klev+1) = 0.0
	551	pstab(jl, klev+1) = 0.0
	552	pstab(jl, 1) = 0.0
	553	pri(jl, klev+1) = 9999.0
	554	ppsi(jl, klev+1) = 0.0
	555	pri(jl, 1) = 0.0
	556	pvph(jl, 1) = 0.0
	557	pulow(jl) = 0.0
	558	pvlow(jl) = 0.0
	559	zulow(jl) = 0.0
	560	zvlow(jl) = 0.0
	561	kkcrith(jl) = klev
	562	kkenvh(jl) = klev
	563	kentp(jl) = klev
	564	kcrit(jl) = 1
	565	ncount(jl) = 0
	566	ll1(jl, klev+1) = .FALSE.
	567	END DO
	568
	569	! * define low-level flow
	570	! ---------------------
	571
	572	DO jk = klev, 2, -1
	573	DO jl = kidia, kfdia
	574	IF (ktest(jl)==1) THEN
	575	zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
	576	prho(jl, jk) = 2.paphm1(jl, jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
	577	pstab(jl, jk) = 2.zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1)) &
	578	(1.-rcpdprho(jl,jk)(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
	579	pstab(jl, jk) = max(pstab(jl,jk), gssec)
	580	END IF
	581	END DO
	582	END DO
	583
	584	! ********************************************************************
	585
	586	! * define blocked flow
	587	! -------------------
	588	DO jk = klev, ilevh, -1
	589	DO jl = kidia, kfdia
	590	IF (jk>=kknub(jl) .AND. jk<=kknul(jl)) THEN
	591	pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	592	pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	593	END IF
	594	END DO
	595	END DO
	596	DO jl = kidia, kfdia
	597	pulow(jl) = pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
	598	pvlow(jl) = pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
	599	znorm(jl) = max(sqrt(pulow(jl)2+pvlow(jl)2), gvsec)
	600	pvph(jl, klev+1) = znorm(jl)
	601	END DO
	602
	603	! ***** setup orography axes and define plane of profiles *****
	604
	605	DO jl = kidia, kfdia
	606	lo = (pulow(jl)<gvsec) .AND. (pulow(jl)>=-gvsec)
	607	IF (lo) THEN
	608	zu = pulow(jl) + 2.*gvsec
	609	ELSE
	610	zu = pulow(jl)
	611	END IF
	612	zphi = atan(pvlow(jl)/zu)
	613	ppsi(jl, klev+1) = ptheta(jl)*rpi/180. - zphi
	614	zb(jl) = 1. - 0.18pgamma(jl) - 0.04pgamma(jl)**2
	615	zc(jl) = 0.48pgamma(jl) + 0.3pgamma(jl)**2
	616	pd1(jl) = zb(jl) - (zb(jl)-zc(jl))(sin(ppsi(jl,klev+1))*2)
	617	pd2(jl) = (zb(jl)-zc(jl))sin(ppsi(jl,klev+1))cos(ppsi(jl,klev+1))
	618	pdmod(jl) = sqrt(pd1(jl)2+pd2(jl)2)
	619	END DO
	620
	621	! ********** define flow in plane of lowlevel stress ***********
	622
	623	DO jk = 1, klev
	624	DO jl = kidia, kfdia
	625	IF (ktest(jl)==1) THEN
	626	zvt1 = pulow(jl)pum1(jl, jk) + pvlow(jl)pvm1(jl, jk)
	627	zvt2 = -pvlow(jl)pum1(jl, jk) + pulow(jl)pvm1(jl, jk)
	628	zvpf(jl, jk) = (zvt1pd1(jl)+zvt2pd2(jl))/(znorm(jl)*pdmod(jl))
	629	END IF
	630	ptau(jl, jk) = 0.0
	631	pzdep(jl, jk) = 0.0
	632	ppsi(jl, jk) = 0.0
	633	ll1(jl, jk) = .FALSE.
	634	END DO
	635	END DO
	636	DO jk = 2, klev
	637	DO jl = kidia, kfdia
	638	IF (ktest(jl)==1) THEN
	639	zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
	640	pvph(jl, jk) = ((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+(papm1(jl, &
	641	jk)-paphm1(jl,jk))*zvpf(jl,jk-1))/zdp(jl, jk)
	642	IF (pvph(jl,jk)<gvsec) THEN
	643	pvph(jl, jk) = gvsec
	644	kcrit(jl) = jk
	645	END IF
	646	END IF
	647	END DO
	648	END DO
	649
	650	! * 2.2 brunt-vaisala frequency and density at half levels.
	651
	652
	653	DO jk = ilevh, klev
	654	DO jl = kidia, kfdia
	655	IF (ktest(jl)==1) THEN
	656	IF (jk>=(kknub(jl)+1) .AND. jk<=kknul(jl)) THEN
	657	zst = zcons2/ptm1(jl, jk)(1.-rcpdprho(jl,jk)*(ptm1(jl, &
	658	jk)-ptm1(jl,jk-1))/zdp(jl,jk))
	659	pstab(jl, klev+1) = pstab(jl, klev+1) + zst*zdp(jl, jk)
	660	pstab(jl, klev+1) = max(pstab(jl,klev+1), gssec)
	661	prho(jl, klev+1) = prho(jl, klev+1) + paphm1(jl, jk)2.zdp(jl, jk) &
	662	*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
	663	END IF
	664	END IF
	665	END DO
	666	END DO
	667
	668	DO jl = kidia, kfdia
	669	pstab(jl, klev+1) = pstab(jl, klev+1)/(papm1(jl,kknul(jl))-papm1(jl,kknub &
	670	(jl)))
	671	prho(jl, klev+1) = prho(jl, klev+1)/(papm1(jl,kknul(jl))-papm1(jl,kknub( &
	672	jl)))
	673	zvar = pstd(jl)
	674	END DO
	675
	676	! * 2.3 mean flow richardson number.
	677	! * and critical height for froude layer
	678
	679
	680	DO jk = 2, klev
	681	DO jl = kidia, kfdia
	682	IF (ktest(jl)==1) THEN
	683	zdwind = max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)), gvsec)
	684	pri(jl, jk) = pstab(jl, jk)(zdp(jl,jk)/(rgprho(jl,jk)zdwind))*2
	685	pri(jl, jk) = max(pri(jl,jk), grcrit)
	686	END IF
	687	END DO
	688	END DO
	689
	690
	691
	692	! * define top of 'envelope' layer
	693	! ----------------------------
	694
	695	DO jl = kidia, kfdia
	696	pnu(jl) = 0.0
	697	znum(jl) = 0.0
	698	END DO
	699
	700	DO jk = 2, klev - 1
	701	DO jl = kidia, kfdia
	702
	703	IF (ktest(jl)==1) THEN
	704
	705	IF (jk>=kknub(jl)) THEN
	706
	707	znum(jl) = pnu(jl)
	708	zwind = (pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/ &
	709	max(sqrt(pulow(jl)2+pvlow(jl)2), gvsec)
	710	zwind = max(sqrt(zwind**2), gvsec)
	711	zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
	712	zstabm = sqrt(max(pstab(jl,jk),gssec))
	713	zstabp = sqrt(max(pstab(jl,jk+1),gssec))
	714	zrhom = prho(jl, jk)
	715	zrhop = prho(jl, jk+1)
	716	pnu(jl) = pnu(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
	717	zwind
	718	IF ((znum(jl)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh( &
	719	jl)==klev)) kkenvh(jl) = jk
	720
	721	END IF
	722
	723	END IF
	724
	725	END DO
	726	END DO
	727
	728	! calculation of a dynamical mixing height for the breaking
	729	! of gravity waves:
	730
	731
	732	DO jl = kidia, kfdia
	733	znup(jl) = 0.0
	734	znum(jl) = 0.0
	735	END DO
	736
	737	DO jk = klev - 1, 2, -1
	738	DO jl = kidia, kfdia
	739
	740	IF (ktest(jl)==1) THEN
	741
	742	znum(jl) = znup(jl)
	743	zwind = (pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/ &
	744	max(sqrt(pulow(jl)2+pvlow(jl)2), gvsec)
	745	zwind = max(sqrt(zwind**2), gvsec)
	746	zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
	747	zstabm = sqrt(max(pstab(jl,jk),gssec))
	748	zstabp = sqrt(max(pstab(jl,jk+1),gssec))
	749	zrhom = prho(jl, jk)
	750	zrhop = prho(jl, jk+1)
	751	znup(jl) = znup(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
	752	zwind
	753	IF ((znum(jl)<=rpi/2.) .AND. (znup(jl)>rpi/2.) .AND. (kkcrith( &
	754	jl)==klev)) kkcrith(jl) = jk
	755
	756	END IF
	757
	758	END DO
	759	END DO
	760
	761	DO jl = kidia, kfdia
	762	kkcrith(jl) = min0(kkcrith(jl), kknu2(jl))
	763	kkcrith(jl) = max0(kkcrith(jl), ilevh*2)
	764	END DO
	765
	766	! directional info for flow blocking *************************
	767
	768	DO jk = ilevh, klev
	769	DO jl = kidia, kfdia
	770	IF (jk>=kkenvh(jl)) THEN
	771	lo = (pum1(jl,jk)<gvsec) .AND. (pum1(jl,jk)>=-gvsec)
	772	IF (lo) THEN
	773	zu = pum1(jl, jk) + 2.*gvsec
	774	ELSE
	775	zu = pum1(jl, jk)
	776	END IF
	777	zphi = atan(pvm1(jl,jk)/zu)
	778	ppsi(jl, jk) = ptheta(jl)*rpi/180. - zphi
	779	END IF
	780	END DO
	781	END DO
	782	! forms the vertical 'leakiness' **************************
	783
	784	alpha = 3.
	785
	786	DO jk = ilevh, klev
	787	DO jl = kidia, kfdia
	788	IF (jk>=kkenvh(jl)) THEN
	789	zggeenv = amax1(1., (pgeom1(jl,kkenvh(jl))+pgeom1(jl, &
	790	kkenvh(jl)-1))/2.)
	791	zggeom1 = amax1(pgeom1(jl,jk), 1.)
	792	zgvar = amax1(pstd(jl)*rg, 1.)
	793	! mod pzdep(jl,jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar))
	794	pzdep(jl, jk) = (pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,jk))/ &
	795	(pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,klev))
	796	END IF
	797	END DO
	798	END DO
	799
	800	RETURN
	801	END SUBROUTINE orosetup
	802	SUBROUTINE gwstress(nlon, nlev, ktest, kcrit, kkenvh, kknu, prho, pstab, &
	803	pvph, pstd, psig, pmea, ppic, ptau, pgeom1, pdmod)
	804
	805	! **** gwstress
	806
	807	! purpose.
	808	! --------
	809
	810	! ** interface.
	811	! ----------
	812	! call gwstress from gwdrag
	813
	814	! explicit arguments :
	815	! --------------------
	816	! ==== inputs ===
	817	! ==== outputs ===
	818
	819	! implicit arguments : none
	820	! --------------------
	821
	822	! method.
	823	! -------
	824
	825
	826	! externals.
	827	! ----------
	828
	829
	830	! reference.
	831	! ----------
	832
	833	! see ecmwf research department documentation of the "i.f.s."
	834
	835	! author.
	836	! -------
	837
	838	! modifications.
	839	! --------------
	840	! f. lott put the new gwd on ifs 22/11/93
	841
	842	! -----------------------------------------------------------------------
	843	USE dimphy
	844	IMPLICIT NONE
	845	include "YOMCST.h"
	846	include "YOEGWD.h"
	847
	848	! -----------------------------------------------------------------------
	849
	850	! * 0.1 arguments
	851	! ---------
	852
	853	INTEGER nlon, nlev
	854	INTEGER kcrit(nlon), ktest(nlon), kkenvh(nlon), kknu(nlon)
	855
	856	REAL prho(nlon, nlev+1), pstab(nlon, nlev+1), ptau(nlon, nlev+1), &
	857	pvph(nlon, nlev+1), pgeom1(nlon, nlev), pstd(nlon)
	858
	859	REAL psig(nlon)
	860	REAL pmea(nlon), ppic(nlon)
	861	REAL pdmod(nlon)
	862
	863	! -----------------------------------------------------------------------
	864
	865	! * 0.2 local arrays
	866	! ------------
	867	INTEGER jl
	868	REAL zblock, zvar, zeff
	869	LOGICAL lo
	870
	871	! -----------------------------------------------------------------------
	872
	873	! * 0.3 functions
	874	! ---------
	875	! ------------------------------------------------------------------
	876
	877	! * 1. initialization
	878	! --------------
	879
	880
	881	! * 3.1 gravity wave stress.
	882
	883
	884
	885	DO jl = kidia, kfdia
	886	IF (ktest(jl)==1) THEN
	887
	888	! effective mountain height above the blocked flow
	889
	890	IF (kkenvh(jl)==klev) THEN
	891	zblock = 0.0
[524]	892	ELSE
[1992]	893	zblock = (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg
	894	END IF
[524]	895
[1992]	896	zvar = ppic(jl) - pmea(jl)
	897	zeff = amax1(0., zvar-zblock)
	898
	899	ptau(jl, klev+1) = prho(jl, klev+1)gkdragpsig(jl)zeff*2/4./ &
	900	pstd(jl)pvph(jl, klev+1)pdmod(jl)*sqrt(pstab(jl,klev+1))
	901
	902	! too small value of stress or low level flow include critical level
	903	! or low level flow: gravity wave stress nul.
	904
	905	lo = (ptau(jl,klev+1)<gtsec) .OR. (kcrit(jl)>=kknu(jl)) .OR. &
	906	(pvph(jl,klev+1)<gvcrit)
	907	! if(lo) ptau(jl,klev+1)=0.0
	908
	909	ELSE
	910
	911	ptau(jl, klev+1) = 0.0
	912
	913	END IF
	914
	915	END DO
	916
	917	RETURN
	918	END SUBROUTINE gwstress
	919	SUBROUTINE gwprofil(nlon, nlev, kgwd, kdx, ktest, kkcrith, kcrit, paphm1, &
	920	prho, pstab, pvph, pri, ptau, pdmod, psig, pvar)
	921
	922	! **** GWPROFIL
	923
	924	! PURPOSE.
	925	! --------
	926
	927	! ** INTERFACE.
	928	! ----------
	929	! FROM GWDRAG
	930
	931	! EXPLICIT ARGUMENTS :
	932	! --------------------
	933	! ==== INPUTS ===
	934	! ==== OUTPUTS ===
	935
	936	! IMPLICIT ARGUMENTS : NONE
	937	! --------------------
	938
	939	! METHOD:
	940	! -------
	941	! THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
	942	! IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
	943	! AND THE TOP OF THE BLOCKED LAYER (KKENVH).
	944	! IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE
	945	! BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR
	946	! NONLINEAR GRAVITY WAVE BREAKING.
	947	! ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
	948	! LEVEL (KCRIT) OR WHEN IT BREAKS.
	949
	950
	951
	952	! EXTERNALS.
	953	! ----------
	954
	955
	956	! REFERENCE.
	957	! ----------
	958
	959	! SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
	960
	961	! AUTHOR.
	962	! -------
	963
	964	! MODIFICATIONS.
	965	! --------------
	966	! PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
	967	! -----------------------------------------------------------------------
	968	USE dimphy
	969	IMPLICIT NONE
	970
	971
	972
	973
	974	include "YOMCST.h"
	975	include "YOEGWD.h"
	976
	977	! -----------------------------------------------------------------------
	978
	979	! * 0.1 ARGUMENTS
	980	! ---------
	981
	982	INTEGER nlon, nlev
	983	INTEGER kkcrith(nlon), kcrit(nlon), kdx(nlon), ktest(nlon)
	984
	985
	986	REAL paphm1(nlon, nlev+1), pstab(nlon, nlev+1), prho(nlon, nlev+1), &
	987	pvph(nlon, nlev+1), pri(nlon, nlev+1), ptau(nlon, nlev+1)
	988
	989	REAL pdmod(nlon), psig(nlon), pvar(nlon)
	990
	991	! -----------------------------------------------------------------------
	992
	993	! * 0.2 LOCAL ARRAYS
	994	! ------------
	995
	996	INTEGER ilevh, ji, kgwd, jl, jk
	997	REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n
	998	REAL zdelp, zdelpt
	999	REAL zdz2(klon, klev), znorm(klon), zoro(klon)
	1000	REAL ztau(klon, klev+1)
	1001
	1002	! -----------------------------------------------------------------------
	1003
	1004	! * 1. INITIALIZATION
	1005	! --------------
	1006
	1007	! print *,' entree gwprofil'
	1008
	1009
	1010	! * COMPUTATIONAL CONSTANTS.
	1011	! ------------- ----------
	1012
	1013	ilevh = klev/3
	1014
	1015	! DO 400 ji=1,kgwd
	1016	! jl=kdx(ji)
	1017	! Modif vectorisation 02/04/2004
	1018	DO jl = kidia, kfdia
	1019	IF (ktest(jl)==1) THEN
	1020	zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl), 1.0)
	1021	ztau(jl, klev+1) = ptau(jl, klev+1)
	1022	END IF
	1023	END DO
	1024
	1025
	1026	DO jk = klev, 2, -1
	1027
	1028	! * 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
	1029	! BLOCKING LAYER.
	1030
	1031	! DO 411 ji=1,kgwd
	1032	! jl=kdx(ji)
	1033	! Modif vectorisation 02/04/2004
	1034	DO jl = kidia, kfdia
	1035	IF (ktest(jl)==1) THEN
	1036	IF (jk>kkcrith(jl)) THEN
	1037	ptau(jl, jk) = ztau(jl, klev+1)
	1038	! ENDIF
	1039	! IF(JK.EQ.KKCRITH(JL)) THEN
[524]	1040	ELSE
[1992]	1041	ptau(jl, jk) = grahilo*ztau(jl, klev+1)
	1042	END IF
	1043	END IF
	1044	END DO
[524]	1045
[1992]	1046	! * 4.15 CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
	1047	! LOW LEVEL FLOW LAYER.
[524]	1048
[1992]	1049
	1050	! * 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
	1051
	1052
	1053	! DO 421 ji=1,kgwd
	1054	! jl=kdx(ji)
	1055	! Modif vectorisation 02/04/2004
	1056	DO jl = kidia, kfdia
	1057	IF (ktest(jl)==1) THEN
	1058	IF (jk<kkcrith(jl)) THEN
	1059	znorm(jl) = gkdragprho(jl, jk)sqrt(pstab(jl,jk))pvph(jl, jk) &
	1060	zoro(jl)
	1061	zdz2(jl, jk) = ptau(jl, jk+1)/max(znorm(jl), gssec)
	1062	END IF
	1063	END IF
	1064	END DO
	1065
	1066	! * 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
	1067	! * AND STRESS: BREAKING EVALUATION AND CRITICAL
	1068	! LEVEL
	1069
	1070
	1071	! DO 431 ji=1,kgwd
	1072	! jl=Kdx(ji)
	1073	! Modif vectorisation 02/04/2004
	1074	DO jl = kidia, kfdia
	1075	IF (ktest(jl)==1) THEN
	1076
	1077	IF (jk<kkcrith(jl)) THEN
	1078	IF ((ptau(jl,jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
	1079	ptau(jl, jk) = 0.0
	1080	ELSE
	1081	zsqr = sqrt(pri(jl,jk))
	1082	zalfa = sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl, jk)
	1083	zriw = pri(jl, jk)(1.-zalfa)/(1+zalfazsqr)**2
	1084	IF (zriw<grcrit) THEN
	1085	zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit
	1086	zb = 1./grcrit + 2./zsqr
	1087	zalpha = 0.5*(-zb+sqrt(zdel))
	1088	zdz2n = (pvph(jl,jk)zalpha)*2/pstab(jl, jk)
	1089	ptau(jl, jk) = znorm(jl)*zdz2n
	1090	ELSE
	1091	ptau(jl, jk) = znorm(jl)*zdz2(jl, jk)
	1092	END IF
	1093	ptau(jl, jk) = min(ptau(jl,jk), ptau(jl,jk+1))
	1094	END IF
	1095	END IF
	1096	END IF
	1097	END DO
	1098
	1099	END DO
	1100
	1101	! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
	1102
	1103	! DO 530 ji=1,kgwd
	1104	! jl=kdx(ji)
	1105	! Modif vectorisation 02/04/2004
	1106	DO jl = kidia, kfdia
	1107	IF (ktest(jl)==1) THEN
	1108	ztau(jl, kkcrith(jl)) = ptau(jl, kkcrith(jl))
	1109	ztau(jl, nstra) = ptau(jl, nstra)
	1110	END IF
	1111	END DO
	1112
	1113	DO jk = 1, klev
	1114
	1115	! DO 532 ji=1,kgwd
	1116	! jl=kdx(ji)
	1117	! Modif vectorisation 02/04/2004
	1118	DO jl = kidia, kfdia
	1119	IF (ktest(jl)==1) THEN
	1120
	1121
	1122	IF (jk>kkcrith(jl)) THEN
	1123
	1124	zdelp = paphm1(jl, jk) - paphm1(jl, klev+1)
	1125	zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1)
	1126	ptau(jl, jk) = ztau(jl, klev+1) + (ztau(jl,kkcrith(jl))-ztau(jl, &
	1127	klev+1))*zdelp/zdelpt
	1128
	1129	END IF
	1130
	1131	END IF
	1132	END DO
	1133
	1134	! REORGANISATION IN THE STRATOSPHERE
	1135
	1136	! DO 533 ji=1,kgwd
	1137	! jl=kdx(ji)
	1138	! Modif vectorisation 02/04/2004
	1139	DO jl = kidia, kfdia
	1140	IF (ktest(jl)==1) THEN
	1141
	1142
	1143	IF (jk<nstra) THEN
	1144
	1145	zdelp = paphm1(jl, nstra)
	1146	zdelpt = paphm1(jl, jk)
	1147	ptau(jl, jk) = ztau(jl, nstra)*zdelpt/zdelp
	1148
	1149	END IF
	1150
	1151	END IF
	1152	END DO
	1153
	1154	! REORGANISATION IN THE TROPOSPHERE
	1155
	1156	! DO 534 ji=1,kgwd
	1157	! jl=kdx(ji)
	1158	! Modif vectorisation 02/04/2004
	1159	DO jl = kidia, kfdia
	1160	IF (ktest(jl)==1) THEN
	1161
	1162
	1163	IF (jk<kkcrith(jl) .AND. jk>nstra) THEN
	1164
	1165	zdelp = paphm1(jl, jk) - paphm1(jl, kkcrith(jl))
	1166	zdelpt = paphm1(jl, nstra) - paphm1(jl, kkcrith(jl))
	1167	ptau(jl, jk) = ztau(jl, kkcrith(jl)) + (ztau(jl,nstra)-ztau(jl, &
	1168	kkcrith(jl)))*zdelp/zdelpt
	1169
	1170	END IF
	1171	END IF
	1172	END DO
	1173
	1174
	1175	END DO
	1176
	1177
	1178	RETURN
	1179	END SUBROUTINE gwprofil
	1180	SUBROUTINE lift_noro(nlon, nlev, dtime, paprs, pplay, plat, pmea, pstd, ppic, &
	1181	ktest, t, u, v, pulow, pvlow, pustr, pvstr, d_t, d_u, d_v)
	1182
	1183	USE dimphy
	1184	IMPLICIT NONE
	1185	! ======================================================================
	1186	! Auteur(s): F.Lott (LMD/CNRS) date: 19950201
	1187	! Objet: Frottement de la montagne Interface
	1188	! ======================================================================
	1189	! Arguments:
	1190	! dtime---input-R- pas d'integration (s)
	1191	! paprs---input-R-pression pour chaque inter-couche (en Pa)
	1192	! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	1193	! t-------input-R-temperature (K)
	1194	! u-------input-R-vitesse horizontale (m/s)
	1195	! v-------input-R-vitesse horizontale (m/s)
	1196
	1197	! d_t-----output-R-increment de la temperature
	1198	! d_u-----output-R-increment de la vitesse u
	1199	! d_v-----output-R-increment de la vitesse v
	1200	! ======================================================================
	1201	include "YOMCST.h"
	1202
	1203	! ARGUMENTS
	1204
	1205	INTEGER nlon, nlev
	1206	REAL dtime
	1207	REAL paprs(klon, klev+1)
	1208	REAL pplay(klon, klev)
	1209	REAL plat(nlon), pmea(nlon)
	1210	REAL pstd(nlon)
	1211	REAL ppic(nlon)
	1212	REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon)
	1213	REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev)
	1214	REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev)
	1215
	1216	INTEGER i, k, ktest(nlon)
	1217
	1218	! Variables locales:
	1219
	1220	REAL zgeom(klon, klev)
	1221	REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev)
	1222	REAL pt(klon, klev), pu(klon, klev), pv(klon, klev)
	1223	REAL papmf(klon, klev), papmh(klon, klev+1)
	1224
	1225	! initialiser les variables de sortie (pour securite)
	1226
	1227	DO i = 1, klon
	1228	pulow(i) = 0.0
	1229	pvlow(i) = 0.0
	1230	pustr(i) = 0.0
	1231	pvstr(i) = 0.0
	1232	END DO
	1233	DO k = 1, klev
	1234	DO i = 1, klon
	1235	d_t(i, k) = 0.0
	1236	d_u(i, k) = 0.0
	1237	d_v(i, k) = 0.0
	1238	pdudt(i, k) = 0.0
	1239	pdvdt(i, k) = 0.0
	1240	pdtdt(i, k) = 0.0
	1241	END DO
	1242	END DO
	1243
	1244	! preparer les variables d'entree (attention: l'ordre des niveaux
	1245	! verticaux augmente du haut vers le bas)
	1246
	1247	DO k = 1, klev
	1248	DO i = 1, klon
	1249	pt(i, k) = t(i, klev-k+1)
	1250	pu(i, k) = u(i, klev-k+1)
	1251	pv(i, k) = v(i, klev-k+1)
	1252	papmf(i, k) = pplay(i, klev-k+1)
	1253	END DO
	1254	END DO
	1255	DO k = 1, klev + 1
	1256	DO i = 1, klon
	1257	papmh(i, k) = paprs(i, klev-k+2)
	1258	END DO
	1259	END DO
	1260	DO i = 1, klon
	1261	zgeom(i, klev) = rdpt(i, klev)log(papmh(i,klev+1)/papmf(i,klev))
	1262	END DO
	1263	DO k = klev - 1, 1, -1
	1264	DO i = 1, klon
	1265	zgeom(i, k) = zgeom(i, k+1) + rd(pt(i,k)+pt(i,k+1))/2.0log(papmf(i,k+ &
	1266	1)/papmf(i,k))
	1267	END DO
	1268	END DO
	1269
	1270	! appeler la routine principale
	1271
	1272	CALL orolift(klon, klev, ktest, dtime, papmh, zgeom, pt, pu, pv, plat, &
	1273	pmea, pstd, ppic, pulow, pvlow, pdudt, pdvdt, pdtdt)
	1274
	1275	DO k = 1, klev
	1276	DO i = 1, klon
	1277	d_u(i, klev+1-k) = dtime*pdudt(i, k)
	1278	d_v(i, klev+1-k) = dtime*pdvdt(i, k)
	1279	d_t(i, klev+1-k) = dtime*pdtdt(i, k)
	1280	pustr(i) = pustr(i) & ! IM BUG .
	1281	! +RGpdudt(i,k)(papmh(i,k+1)-papmh(i,k))
	1282	+pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
	1283	pvstr(i) = pvstr(i) & ! IM BUG .
	1284	! +RGpdvdt(i,k)(papmh(i,k+1)-papmh(i,k))
	1285	+pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
	1286	END DO
	1287	END DO
	1288
	1289	RETURN
	1290	END SUBROUTINE lift_noro
	1291	SUBROUTINE orolift(nlon, nlev, ktest, ptsphy, paphm1, pgeom1, ptm1, pum1, &
	1292	pvm1, plat, pmea, pvaror, ppic & ! OUTPUTS
	1293	, pulow, pvlow, pvom, pvol, pte)
	1294
	1295
	1296	! **** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
	1297
	1298	! PURPOSE.
	1299	! --------
	1300
	1301	! ** INTERFACE.
	1302	! ----------
	1303	! CALLED FROM *lift_noro
	1304	! ----------
	1305
	1306	! AUTHOR.
	1307	! -------
	1308	! F.LOTT LMD 22/11/95
	1309
	1310	USE dimphy
	1311	IMPLICIT NONE
	1312
	1313
	1314	include "YOMCST.h"
	1315	include "YOEGWD.h"
	1316	! -----------------------------------------------------------------------
	1317
	1318	! * 0.1 ARGUMENTS
	1319	! ---------
	1320
	1321
	1322	INTEGER nlon, nlev
	1323	REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(nlon), &
	1324	pvlow(nlon)
	1325	REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), plat(nlon), &
	1326	pmea(nlon), pvaror(nlon), ppic(nlon), pgeom1(nlon, nlev), &
	1327	paphm1(nlon, nlev+1)
	1328
	1329	INTEGER ktest(nlon)
	1330	REAL ptsphy
	1331	! -----------------------------------------------------------------------
	1332
	1333	! * 0.2 LOCAL ARRAYS
	1334	! ------------
	1335	LOGICAL lifthigh
	1336	! ym integer klevm1, jl, ilevh, jk
	1337	INTEGER jl, ilevh, jk
	1338	REAL zcons1, ztmst, zrtmst, zpi, zhgeo
	1339	REAL zdelp, zslow, zsqua, zscav, zbet
	1340	INTEGER iknub(klon), iknul(klon)
	1341	LOGICAL ll1(klon, klev+1)
	1342
	1343	REAL ztau(klon, klev+1), ztav(klon, klev+1), zrho(klon, klev+1)
	1344	REAL zdudt(klon), zdvdt(klon)
	1345	REAL zhcrit(klon, klev)
	1346	CHARACTER (LEN=20) :: modname = 'orografi'
	1347	CHARACTER (LEN=80) :: abort_message
	1348	! -----------------------------------------------------------------------
	1349
	1350	! * 1.1 INITIALIZATIONS
	1351	! ---------------
	1352
	1353	lifthigh = .FALSE.
	1354
	1355	IF (nlon/=klon .OR. nlev/=klev) THEN
	1356	abort_message = 'pb dimension'
[2408]	1357	CALL abort_physic(modname, abort_message, 1)
[1992]	1358	END IF
	1359	zcons1 = 1./rd
	1360	! ym KLEVM1=KLEV-1
	1361	ztmst = ptsphy
	1362	zrtmst = 1./ztmst
	1363	zpi = acos(-1.)
	1364
	1365	DO jl = kidia, kfdia
	1366	zrho(jl, klev+1) = 0.0
	1367	pulow(jl) = 0.0
	1368	pvlow(jl) = 0.0
	1369	iknub(jl) = klev
	1370	iknul(jl) = klev
	1371	ilevh = klev/3
	1372	ll1(jl, klev+1) = .FALSE.
	1373	DO jk = 1, klev
	1374	pvom(jl, jk) = 0.0
	1375	pvol(jl, jk) = 0.0
	1376	pte(jl, jk) = 0.0
	1377	END DO
	1378	END DO
	1379
	1380
	1381	! * 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
	1382	! * LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
	1383	! * THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
	1384
	1385
	1386
	1387	DO jk = klev, 1, -1
	1388	DO jl = kidia, kfdia
	1389	IF (ktest(jl)==1) THEN
	1390	zhcrit(jl, jk) = amax1(ppic(jl)-pmea(jl), 100.)
	1391	zhgeo = pgeom1(jl, jk)/rg
	1392	ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
	1393	IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
	1394	iknub(jl) = jk
	1395	END IF
	1396	END IF
	1397	END DO
	1398	END DO
	1399
	1400	DO jl = kidia, kfdia
	1401	IF (ktest(jl)==1) THEN
	1402	iknub(jl) = max(iknub(jl), klev/2)
	1403	iknul(jl) = max(iknul(jl), 2*klev/3)
	1404	IF (iknub(jl)>nktopg) iknub(jl) = nktopg
	1405	IF (iknub(jl)==nktopg) iknul(jl) = klev
	1406	IF (iknub(jl)==iknul(jl)) iknub(jl) = iknul(jl) - 1
	1407	END IF
	1408	END DO
	1409
	1410	! do 2011 jl=kidia,kfdia
	1411	! IF(KTEST(JL).EQ.1) THEN
	1412	! print *,' iknul= ',iknul(jl),' iknub=',iknub(jl)
	1413	! ENDIF
	1414	! 2011 continue
	1415
	1416	! PRINT *,' DANS OROLIFT: 2010'
	1417
	1418	DO jk = klev, 2, -1
	1419	DO jl = kidia, kfdia
	1420	zrho(jl, jk) = 2.paphm1(jl, jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
	1421	END DO
	1422	END DO
	1423	! PRINT *,' DANS OROLIFT: 223'
	1424
	1425	! ********************************************************************
	1426
	1427	! * DEFINE LOW LEVEL FLOW
	1428	! -------------------
	1429	DO jk = klev, 1, -1
	1430	DO jl = kidia, kfdia
	1431	IF (ktest(jl)==1) THEN
	1432	IF (jk>=iknub(jl) .AND. jk<=iknul(jl)) THEN
	1433	pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) &
	1434	)
	1435	pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) &
	1436	)
	1437	zrho(jl, klev+1) = zrho(jl, klev+1) + zrho(jl, jk)*(paphm1(jl,jk+1) &
	1438	-paphm1(jl,jk))
	1439	END IF
	1440	END IF
	1441	END DO
	1442	END DO
	1443	DO jl = kidia, kfdia
	1444	IF (ktest(jl)==1) THEN
	1445	pulow(jl) = pulow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
	1446	pvlow(jl) = pvlow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
	1447	zrho(jl, klev+1) = zrho(jl, klev+1)/(paphm1(jl,iknul(jl)+1)-paphm1(jl, &
	1448	iknub(jl)))
	1449	END IF
	1450	END DO
	1451
	1452	! ***********************************************************
	1453
	1454	! * 3. COMPUTE MOUNTAIN LIFT
	1455
	1456
	1457	DO jl = kidia, kfdia
	1458	IF (ktest(jl)==1) THEN
	1459	ztau(jl, klev+1) = -gkliftzrho(jl, klev+1)2.romega & ! *
	1460	! (2PVAROR(JL)+PMEA(JL))
	1461	2pvaror(jl)sin(zpi/180.plat(jl))pvlow(jl)
	1462	ztav(jl, klev+1) = gkliftzrho(jl, klev+1)2.romega & ! *
	1463	! (2PVAROR(JL)+PMEA(JL))
	1464	2pvaror(jl)sin(zpi/180.plat(jl))pulow(jl)
	1465	ELSE
	1466	ztau(jl, klev+1) = 0.0
	1467	ztav(jl, klev+1) = 0.0
	1468	END IF
	1469	END DO
	1470
	1471	! * 4. COMPUTE LIFT PROFILE
	1472	! * --------------------
	1473
	1474
	1475
	1476	DO jk = 1, klev
	1477	DO jl = kidia, kfdia
	1478	IF (ktest(jl)==1) THEN
	1479	ztau(jl, jk) = ztau(jl, klev+1)*paphm1(jl, jk)/paphm1(jl, klev+1)
	1480	ztav(jl, jk) = ztav(jl, klev+1)*paphm1(jl, jk)/paphm1(jl, klev+1)
[524]	1481	ELSE
[1992]	1482	ztau(jl, jk) = 0.0
	1483	ztav(jl, jk) = 0.0
	1484	END IF
	1485	END DO
	1486	END DO
[524]	1487
	1488
[1992]	1489	! * 5. COMPUTE TENDENCIES.
	1490	! * -------------------
	1491	IF (lifthigh) THEN
	1492	! PRINT *,' DANS OROLIFT: 500'
[524]	1493
[1992]	1494	! EXPLICIT SOLUTION AT ALL LEVELS
[766]	1495
[1992]	1496	DO jk = 1, klev
	1497	DO jl = kidia, kfdia
	1498	IF (ktest(jl)==1) THEN
	1499	zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
	1500	zdudt(jl) = -rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
	1501	zdvdt(jl) = -rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
	1502	END IF
	1503	END DO
	1504	END DO
[766]	1505
[1992]	1506	! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
[766]	1507
[1992]	1508	DO jk = 1, klev
	1509	DO jl = kidia, kfdia
	1510	IF (ktest(jl)==1) THEN
[766]	1511
[1992]	1512	zslow = sqrt(pulow(jl)2+pvlow(jl)2)
	1513	zsqua = amax1(sqrt(pum1(jl,jk)2+pvm1(jl,jk)2), gvsec)
	1514	zscav = -zdudt(jl)pvm1(jl, jk) + zdvdt(jl)pum1(jl, jk)
	1515	IF (zsqua>gvsec) THEN
	1516	pvom(jl, jk) = -zscavpvm1(jl, jk)/zsqua*2
	1517	pvol(jl, jk) = zscavpum1(jl, jk)/zsqua*2
	1518	ELSE
	1519	pvom(jl, jk) = 0.0
	1520	pvol(jl, jk) = 0.0
	1521	END IF
	1522	zsqua = sqrt(pum1(jl,jk)2+pum1(jl,jk)2)
	1523	IF (zsqua<zslow) THEN
	1524	pvom(jl, jk) = zsqua/zslow*pvom(jl, jk)
	1525	pvol(jl, jk) = zsqua/zslow*pvol(jl, jk)
	1526	END IF
[766]	1527
[1992]	1528	END IF
	1529	END DO
	1530	END DO
	1531
	1532	! 6. LOW LEVEL LIFT, SEMI IMPLICIT:
	1533	! ----------------------------------
	1534
	1535	ELSE
	1536
	1537	DO jl = kidia, kfdia
	1538	IF (ktest(jl)==1) THEN
	1539	DO jk = klev, iknub(jl), -1
	1540	zbet = gklift2.romegasin(zpi/180.plat(jl))ztmst &
	1541	(pgeom1(jl,iknub(jl)-1)-pgeom1(jl,jk))/ &
	1542	(pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
	1543	zdudt(jl) = -pum1(jl, jk)/ztmst/(1+zbet**2)
	1544	zdvdt(jl) = -pvm1(jl, jk)/ztmst/(1+zbet**2)
	1545	pvom(jl, jk) = zbet*2zdudt(jl) - zbet*zdvdt(jl)
	1546	pvol(jl, jk) = zbetzdudt(jl) + zbet2zdvdt(jl)
	1547	END DO
	1548	END IF
	1549	END DO
	1550
	1551	END IF
	1552
	1553	RETURN
	1554	END SUBROUTINE orolift
	1555
	1556
	1557	SUBROUTINE sugwd(nlon, nlev, paprs, pplay)
	1558	USE dimphy
	1559	USE mod_phys_lmdz_para
	1560	USE mod_grid_phy_lmdz
	1561	! USE parallel
	1562
	1563	! **** SUGWD INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
	1564
	1565	! PURPOSE.
	1566	! --------
	1567	! INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
	1568	! GRAVITY WAVE DRAG PARAMETRIZATION.
	1569
	1570	! ** INTERFACE.
	1571	! ----------
	1572	! CALL SUGWD FROM SUPHEC
	1573	! ----- ------
	1574
	1575	! EXPLICIT ARGUMENTS :
	1576	! --------------------
	1577	! PSIG : VERTICAL COORDINATE TABLE
	1578	! NLEV : NUMBER OF MODEL LEVELS
	1579
	1580	! IMPLICIT ARGUMENTS :
	1581	! --------------------
	1582	! COMMON YOEGWD
	1583
	1584	! METHOD.
	1585	! -------
	1586	! SEE DOCUMENTATION
	1587
	1588	! EXTERNALS.
	1589	! ----------
	1590	! NONE
	1591
	1592	! REFERENCE.
	1593	! ----------
	1594	! ECMWF Research Department documentation of the IFS
	1595
	1596	! AUTHOR.
	1597	! -------
	1598	! MARTIN MILLER ECMWF
	1599
	1600	! MODIFICATIONS.
	1601	! --------------
	1602	! ORIGINAL : 90-01-01
	1603	! ------------------------------------------------------------------
	1604	IMPLICIT NONE
	1605
	1606	! -----------------------------------------------------------------
	1607	include "YOEGWD.h"
	1608	! ----------------------------------------------------------------
	1609
	1610	INTEGER nlon, nlev, jk
	1611	REAL paprs(nlon, nlev+1)
	1612	REAL pplay(nlon, nlev)
	1613	REAL zpr, zstra, zsigt, zpm1r
	1614	REAL :: pplay_glo(klon_glo, nlev)
	1615	REAL :: paprs_glo(klon_glo, nlev+1)
	1616
	1617	! * 1. SET THE VALUES OF THE PARAMETERS
	1618	! --------------------------------
	1619
	1620
	1621	PRINT *, ' DANS SUGWD NLEV=', nlev
	1622	ghmax = 10000.
	1623
	1624	zpr = 100000.
	1625	zstra = 0.1
	1626	zsigt = 0.94
	1627	! old ZPR=80000.
	1628	! old ZSIGT=0.85
	1629
	1630
	1631	CALL gather(pplay, pplay_glo)
	1632	CALL bcast(pplay_glo)
	1633	CALL gather(paprs, paprs_glo)
	1634	CALL bcast(paprs_glo)
	1635
	1636
	1637	DO jk = 1, nlev
	1638	zpm1r = pplay_glo((klon_glo/2)+1, jk)/paprs_glo((klon_glo/2)+1, 1)
	1639	IF (zpm1r>=zsigt) THEN
	1640	nktopg = jk
	1641	END IF
	1642	zpm1r = pplay_glo((klon_glo/2)+1, jk)/paprs_glo((klon_glo/2)+1, 1)
	1643	IF (zpm1r>=zstra) THEN
	1644	nstra = jk
	1645	END IF
	1646	END DO
	1647
	1648
	1649
	1650	! inversion car dans orodrag on compte les niveaux a l'envers
	1651	nktopg = nlev - nktopg + 1
	1652	nstra = nlev - nstra
	1653	PRINT *, ' DANS SUGWD nktopg=', nktopg
	1654	PRINT *, ' DANS SUGWD nstra=', nstra
	1655
	1656	gsigcr = 0.80
	1657
[2408]	1658	! Values now specified in run.def, or conf_phys_m.F90
	1659	! gkdrag = 0.2
	1660	! grahilo = 1.
	1661	! grcrit = 0.01
	1662	! gfrcrit = 1.0
	1663	! gkwake = 0.50
	1664	! gklift = 0.50
[1992]	1665	gvcrit = 0.0
	1666
	1667	! ----------------------------------------------------------------
	1668
	1669	! * 2. SET VALUES OF SECURITY PARAMETERS
	1670	! ---------------------------------
	1671
	1672
	1673	gvsec = 0.10
	1674	gssec = 1.E-12
	1675
	1676	gtsec = 1.E-07
	1677
	1678	! ----------------------------------------------------------------
	1679
	1680	RETURN
	1681	END SUBROUTINE sugwd

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

Download in other formats:

Original Format