Context Navigation

orografi.F @ 281

Last change on this file since 281 was 177, checked in by lmdzadmin, 24 years ago

Lots of stuff, plus particulierement:

appel a ORCHIDEE en etat de marche (pb de grille subsiste)
modifs de Pascale sur soil dans le cas ou ok_veget=false -

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 46.4 KB

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

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

Download in other formats:

Original Format