Context Navigation

orografi.F90 @ 2306

Last change on this file since 2306 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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: 44.0 KB

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

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

Download in other formats:

Original Format