Context Navigation

orografi.F @ 2019

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

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

Download in other formats:

Original Format