Context Navigation

orografi.F @ 4944

Last change on this file since 4944 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 48.8 KB

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

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

Download in other formats:

Original Format