Context Navigation

orografi.F @ 2828

Last change on this file since 2828 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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

File size: 48.5 KB

Rev	Line
[524]	1	!
[1299]	2	! $Id: orografi.F 1299 2010-01-20 14:27:21Z fhourdin $
[524]	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
[766]	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======================================================================
[766]	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
[766]	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
[766]	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
[766]	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
[766]	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-----------------------------------------------------------------------
[766]	471	USE dimphy
[524]	472	implicit none
	473	c
	474
[766]	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-----------------------------------------------------------------------
[766]	940	USE dimphy
[524]	941	implicit none
[766]	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-----------------------------------------------------------------------
[766]	1075	USE dimphy
[524]	1076	implicit none
	1077	C
	1078
	1079	C
	1080
[766]	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
[766]	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======================================================================
[766]	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
[766]	1448	USE dimphy
[524]	1449	implicit none
	1450	C
	1451	C
[766]	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 ------------
[776]	1483	logical lifthigh
[766]	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)
[1299]	1499	CHARACTER (LEN=20) :: modname='orografi'
	1500	CHARACTER (LEN=80) :: abort_message
[524]	1501	C-----------------------------------------------------------------------
	1502	C
	1503	C* 1.1 INITIALIZATIONS
	1504	C ---------------
	1505
	1506	LIFTHIGH=.FALSE.
	1507
[1299]	1508	IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)THEN
	1509	abort_message = 'pb dimension'
	1510	CALL abort_gcm (modname,abort_message,1)
	1511	ENDIF
[524]	1512	ZCONS1=1./RD
[766]	1513	cym KLEVM1=KLEV-1
[524]	1514	ZTMST=PTSPHY
	1515	ZRTMST=1./ZTMST
	1516	ZPI=ACOS(-1.)
	1517	C
	1518	DO 1001 JL=kidia,kfdia
	1519	ZRHO(JL,KLEV+1) =0.0
	1520	PULOW(JL) =0.0
	1521	PVLOW(JL) =0.0
	1522	iknub(JL) =klev
	1523	iknul(JL) =klev
	1524	ilevh=klev/3
	1525	ll1(jl,klev+1)=.false.
	1526	DO 1000 JK=1,KLEV
	1527	PVOM(JL,JK)=0.0
	1528	PVOL(JL,JK)=0.0
	1529	PTE (JL,JK)=0.0
	1530	1000 CONTINUE
	1531	1001 CONTINUE
	1532
	1533	C
	1534	C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
	1535	C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
	1536	C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
	1537	C
	1538	C
	1539	C
	1540	DO 2006 JK=KLEV,1,-1
	1541	DO 2007 JL=kidia,kfdia
	1542	IF(KTEST(JL).EQ.1) THEN
	1543	ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.)
	1544	ZHGEO=PGEOM1(JL,JK)/RG
	1545	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
[675]	1546	IF(ll1(JL,JK).neqv.ll1(JL,JK+1)) THEN
[524]	1547	iknub(JL)=JK
	1548	ENDIF
	1549	ENDIF
	1550	2007 CONTINUE
	1551	2006 CONTINUE
	1552	C
	1553	do 2010 jl=kidia,kfdia
	1554	IF(KTEST(JL).EQ.1) THEN
	1555	iknub(jl)=max(iknub(jl),klev/2)
	1556	iknul(jl)=max(iknul(jl),2*klev/3)
	1557	if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
	1558	if(iknub(jl).eq.nktopg) iknul(jl)=klev
	1559	if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
	1560	ENDIF
	1561	2010 continue
	1562
	1563	C do 2011 jl=kidia,kfdia
	1564	C IF(KTEST(JL).EQ.1) THEN
	1565	C print *,' iknul= ',iknul(jl),' iknub=',iknub(jl)
	1566	C ENDIF
	1567	C2011 continue
	1568
	1569	C PRINT *,' DANS OROLIFT: 2010'
	1570
	1571	DO 223 JK=KLEV,2,-1
	1572	DO 222 JL=kidia,kfdia
	1573	ZRHO(JL,JK)=2.PAPHM1(JL,JK)ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
	1574	222 CONTINUE
	1575	223 CONTINUE
	1576	C PRINT *,' DANS OROLIFT: 223'
	1577
	1578	C********************************************************************
	1579	C
	1580	C* DEFINE LOW LEVEL FLOW
	1581	C -------------------
	1582	DO 2115 JK=klev,1,-1
	1583	DO 2116 JL=kidia,kfdia
	1584	IF(KTEST(JL).EQ.1) THEN
	1585	if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
	1586	pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
	1587	pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
	1588	zrho(JL,klev+1)=zrho(JL,klev+1)
	1589	* +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
	1590	end if
	1591	ENDIF
	1592	2116 CONTINUE
	1593	2115 CONTINUE
	1594	DO 2110 JL=kidia,kfdia
	1595	IF(KTEST(JL).EQ.1) THEN
	1596	pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
	1597	pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
	1598	zrho(JL,klev+1)=zrho(JL,klev+1)
	1599	* /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
	1600	ENDIF
	1601	2110 CONTINUE
	1602
	1603
	1604	200 CONTINUE
	1605
	1606	C***********************************************************
	1607	C
	1608	C* 3. COMPUTE MOUNTAIN LIFT
	1609	C
	1610	300 CONTINUE
	1611	C
	1612	DO 301 JL=kidia,kfdia
	1613	IF(KTEST(JL).EQ.1) THEN
	1614	ZTAU(JL,KLEV+1)= - GKLIFTZRHO(JL,KLEV+1)2.ROMEGA
	1615	C * (2PVAROR(JL)+PMEA(JL))
	1616	* 2PVAROR(JL)
	1617	* SIN(ZPI/180.PLAT(JL))PVLOW(JL)
	1618	ZTAV(JL,KLEV+1)= GKLIFTZRHO(JL,KLEV+1)2.ROMEGA
	1619	C * (2PVAROR(JL)+PMEA(JL))
	1620	* 2PVAROR(JL)
	1621	* SIN(ZPI/180.PLAT(JL))PULOW(JL)
	1622	ELSE
	1623	ZTAU(JL,KLEV+1)=0.0
	1624	ZTAV(JL,KLEV+1)=0.0
	1625	ENDIF
	1626	301 CONTINUE
	1627
	1628	C
	1629	C* 4. COMPUTE LIFT PROFILE
	1630	C* --------------------
	1631	C
	1632
	1633	400 CONTINUE
	1634
	1635	DO 401 JK=1,KLEV
	1636	DO 401 JL=kidia,kfdia
	1637	IF(KTEST(JL).EQ.1) THEN
	1638	ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
	1639	ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
	1640	ELSE
	1641	ZTAU(JL,JK)=0.0
	1642	ZTAV(JL,JK)=0.0
	1643	ENDIF
	1644	401 CONTINUE
	1645	C
	1646	C
	1647	C* 5. COMPUTE TENDENCIES.
	1648	C* -------------------
	1649	IF(LIFTHIGH)THEN
	1650	C
	1651	500 CONTINUE
	1652	C PRINT *,' DANS OROLIFT: 500'
	1653	C
	1654	C EXPLICIT SOLUTION AT ALL LEVELS
	1655	C
	1656	DO 524 JK=1,klev
	1657	DO 523 JL=KIDIA,KFDIA
	1658	IF(KTEST(JL).EQ.1) THEN
	1659	ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
	1660	ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP
	1661	ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP
	1662	ENDIF
	1663	523 CONTINUE
	1664	524 CONTINUE
	1665	C
	1666	C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
	1667	C
	1668	DO 530 JK=1,klev
	1669	DO 530 JL=KIDIA,KFDIA
	1670	IF(KTEST(JL).EQ.1) THEN
	1671
	1672	ZSLOW=SQRT(PULOW(JL)2+PVLOW(JL)2)
	1673	ZSQUA=AMAX1(SQRT(PUM1(JL,JK)2+PVM1(JL,JK)2),GVSEC)
	1674	ZSCAV=-ZDUDT(JL)PVM1(JL,JK)+ZDVDT(JL)PUM1(JL,JK)
	1675	IF(ZSQUA.GT.GVSEC)THEN
	1676	PVOM(JL,JK)=-ZSCAVPVM1(JL,JK)/ZSQUA*2
	1677	PVOL(JL,JK)= ZSCAVPUM1(JL,JK)/ZSQUA*2
	1678	ELSE
	1679	PVOM(JL,JK)=0.0
	1680	PVOL(JL,JK)=0.0
	1681	ENDIF
	1682	ZSQUA=SQRT(PUM1(JL,JK)2+PUM1(JL,JK)2)
	1683	IF(ZSQUA.LT.ZSLOW)THEN
	1684	PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK)
	1685	PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK)
	1686	ENDIF
	1687
	1688	ENDIF
	1689	530 CONTINUE
	1690	C
	1691	C 6. LOW LEVEL LIFT, SEMI IMPLICIT:
	1692	C ----------------------------------
	1693
	1694	ELSE
	1695
	1696	DO 601 JL=KIDIA,KFDIA
	1697	IF(KTEST(JL).EQ.1) THEN
	1698	DO JK=KLEV,IKNUB(JL),-1
	1699	ZBET=GKLIFT2.ROMEGASIN(ZPI/180.PLAT(JL))ztmst
	1700	* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL, JK))/
	1701	* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV))
	1702	ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2)
	1703	ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2)
	1704	PVOM(JL,JK)= ZBET*2ZDUDT(JL) - ZBET *ZDVDT(JL)
	1705	PVOL(JL,JK)= ZBET ZDUDT(JL) + ZBET2ZDVDT(JL)
	1706	ENDDO
	1707	ENDIF
	1708	601 CONTINUE
	1709
	1710	ENDIF
	1711
	1712	RETURN
	1713	END
[766]	1714
	1715
[524]	1716	SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay)
[766]	1717	USE dimphy
[776]	1718	USE mod_phys_lmdz_para
	1719	USE mod_grid_phy_lmdz
[766]	1720	c USE parallel
[524]	1721	C
	1722	C**** SUGWD INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
	1723	C
	1724	C PURPOSE.
	1725	C --------
	1726	C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
	1727	C GRAVITY WAVE DRAG PARAMETRIZATION.
	1728	C
	1729	C** INTERFACE.
	1730	C ----------
	1731	C CALL SUGWD FROM SUPHEC
	1732	C ----- ------
	1733	C
	1734	C EXPLICIT ARGUMENTS :
	1735	C --------------------
	1736	C PSIG : VERTICAL COORDINATE TABLE
	1737	C NLEV : NUMBER OF MODEL LEVELS
	1738	C
	1739	C IMPLICIT ARGUMENTS :
	1740	C --------------------
	1741	C COMMON YOEGWD
	1742	C
	1743	C METHOD.
	1744	C -------
	1745	C SEE DOCUMENTATION
	1746	C
	1747	C EXTERNALS.
	1748	C ----------
	1749	C NONE
	1750	C
	1751	C REFERENCE.
	1752	C ----------
	1753	C ECMWF Research Department documentation of the IFS
	1754	C
	1755	C AUTHOR.
	1756	C -------
	1757	C MARTIN MILLER ECMWF
	1758	C
	1759	C MODIFICATIONS.
	1760	C --------------
	1761	C ORIGINAL : 90-01-01
	1762	C ------------------------------------------------------------------
	1763	implicit none
	1764	C
	1765	C -----------------------------------------------------------------
	1766	#include "YOEGWD.h"
	1767	C ----------------------------------------------------------------
	1768	C
	1769	integer nlon,nlev, jk
	1770	REAL paprs(nlon,nlev+1)
	1771	REAL pplay(nlon,nlev)
	1772	real zpr,zstra,zsigt,zpm1r
[776]	1773	REAL :: pplay_glo(klon_glo,nlev)
	1774	REAL :: paprs_glo(klon_glo,nlev+1)
[766]	1775
[524]	1776	C
	1777	C* 1. SET THE VALUES OF THE PARAMETERS
	1778	C --------------------------------
	1779	C
	1780	100 CONTINUE
	1781	C
	1782	PRINT *,' DANS SUGWD NLEV=',NLEV
	1783	GHMAX=10000.
	1784	C
	1785	ZPR=100000.
	1786	ZSTRA=0.1
	1787	ZSIGT=0.94
	1788	cold ZPR=80000.
	1789	cold ZSIGT=0.85
	1790	C
[766]	1791
[776]	1792	CALL gather(pplay,pplay_glo)
	1793	CALL bcast(pplay_glo)
	1794	CALL gather(paprs,paprs_glo)
	1795	CALL bcast(paprs_glo)
	1796
	1797
[524]	1798	DO 110 JK=1,NLEV
[1279]	1799	ZPM1R=pplay_glo((klon_glo/2)+1,jk)/paprs_glo((klon_glo/2)+1,1)
[524]	1800	IF(ZPM1R.GE.ZSIGT)THEN
	1801	nktopg=JK
	1802	ENDIF
[1279]	1803	ZPM1R=pplay_glo((klon_glo/2)+1,jk)/paprs_glo((klon_glo/2)+1,1)
[524]	1804	IF(ZPM1R.GE.ZSTRA)THEN
	1805	NSTRA=JK
	1806	ENDIF
	1807	110 CONTINUE
[766]	1808
	1809
[524]	1810	c
	1811	c inversion car dans orodrag on compte les niveaux a l'envers
	1812	nktopg=nlev-nktopg+1
	1813	nstra=nlev-nstra
	1814	print *,' DANS SUGWD nktopg=', nktopg
	1815	print *,' DANS SUGWD nstra=', nstra
	1816	C
	1817	GSIGCR=0.80
	1818	C
	1819	GKDRAG=0.2
	1820	GRAHILO=1.
	1821	GRCRIT=0.01
	1822	GFRCRIT=1.0
	1823	GKWAKE=0.50
	1824	C
	1825	GKLIFT=0.50
	1826	GVCRIT =0.0
	1827	C
	1828	C
	1829	C ----------------------------------------------------------------
	1830	C
	1831	C* 2. SET VALUES OF SECURITY PARAMETERS
	1832	C ---------------------------------
	1833	C
	1834	200 CONTINUE
	1835	C
	1836	GVSEC=0.10
	1837	GSSEC=1.E-12
	1838	C
	1839	GTSEC=1.E-07
	1840	C
	1841	C ----------------------------------------------------------------
	1842	C
	1843	RETURN
	1844	END

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

Download in other formats:

Original Format