Context Navigation

orografi_strato.F @ 3702

Last change on this file since 3702 was 1491, checked in by musat, 14 years ago
Couplage Thermiques/Emanuel? seulement si T<0 avec iflag_coupl==2 (FH) Correction bug initialisation ptconvth (FH/JYG) Correction dimensions tableau paprs_glo definition ZPM1R pour le 1D (JYG) IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 57.6 KB

Rev	Line
[1001]	1	SUBROUTINE drag_noro_strato (nlon,nlev,dtime,paprs,pplay,
	2	e pmea,pstd, psig, pgam, pthe,ppic,pval,
	3	e kgwd,kdx,ktest,
	4	e t, u, v,
	5	s pulow, pvlow, pustr, pvstr,
	6	s d_t, d_u, d_v)
	7	c
	8	USE dimphy
	9	IMPLICIT none
	10	c======================================================================
	11	c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
	12	c Object: Mountain drag interface. Made necessary because:
	13	C 1. in the LMD-GCM Layers are from bottom to top,
	14	C contrary to most European GCM.
	15	c 2. the altitude above ground of each model layers
	16	c needs to be known (variable zgeom)
	17	c======================================================================
	18	c Explicit Arguments:
	19	c ==================
	20	c nlon----input-I-Total number of horizontal points that get into physics
	21	c nlev----input-I-Number of vertical levels
	22	c dtime---input-R-Time-step (s)
	23	c paprs---input-R-Pressure in semi layers (Pa)
	24	c pplay---input-R-Pressure model-layers (Pa)
	25	c t-------input-R-temperature (K)
	26	c u-------input-R-Horizontal wind (m/s)
	27	c v-------input-R-Meridional wind (m/s)
	28	c pmea----input-R-Mean Orography (m)
	29	C pstd----input-R-SSO standard deviation (m)
	30	c psig----input-R-SSO slope
	31	c pgam----input-R-SSO Anisotropy
	32	c pthe----input-R-SSO Angle
	33	c ppic----input-R-SSO Peacks elevation (m)
	34	c pval----input-R-SSO Valleys elevation (m)
	35	c
	36	c kgwd- -input-I: Total nb of points where the orography schemes are active
	37	c ktest--input-I: Flags to indicate active points
	38	c kdx----input-I: Locate the physical location of an active point.
	39
	40	c pulow, pvlow -output-R: Low-level wind
	41	c pustr, pvstr -output-R: Surface stress due to SSO drag (Pa)
	42	c
	43	c d_t-----output-R: T increment
	44	c d_u-----output-R: U increment
	45	c d_v-----output-R: V increment
	46	c
	47	c Implicit Arguments:
	48	c ===================
	49	c
	50	c iim--common-I: Number of longitude intervals
	51	c jjm--common-I: Number of latitude intervals
	52	c klon-common-I: Number of points seen by the physics
	53	c (iim+1)*(jjm+1) for instance
	54	c klev-common-I: Number of vertical layers
	55	c======================================================================
	56	c Local Variables:
	57	c ================
	58	c
	59	c zgeom-----R: Altitude of layer above ground
	60	c pt, pu, pv --R: t u v from top to bottom
	61	c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
	62	c papmf: pressure at model layer (from top to bottom)
	63	c papmh: pressure at model 1/2 layer (from top to bottom)
	64	c
	65	c======================================================================
	66	cym#include "dimensions.h"
	67	cym#include "dimphy.h"
	68	#include "YOMCST.h"
	69	#include "YOEGWD.h"
	70	c
	71	c ARGUMENTS
	72	c
	73	INTEGER nlon,nlev
	74	REAL dtime
	75	REAL paprs(nlon,nlev+1)
	76	REAL pplay(nlon,nlev)
	77	REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
	78	REAL ppic(nlon),pval(nlon)
	79	REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
	80	REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
	81	REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
	82	c
	83	INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
	84	c
	85	c LOCAL VARIABLES:
	86	c
	87	REAL zgeom(klon,klev)
	88	REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
	89	REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
	90	REAL papmf(klon,klev),papmh(klon,klev+1)
[1403]	91	CHARACTER (LEN=20) :: modname='orografi_strato'
	92	CHARACTER (LEN=80) :: abort_message
[1001]	93	c
	94	c INITIALIZE OUTPUT VARIABLES
	95	c
	96	DO i = 1,klon
	97	pulow(i) = 0.0
	98	pvlow(i) = 0.0
	99	pustr(i) = 0.0
	100	pvstr(i) = 0.0
	101	ENDDO
	102	DO k = 1, klev
	103	DO i = 1, klon
	104	d_t(i,k) = 0.0
	105	d_u(i,k) = 0.0
	106	d_v(i,k) = 0.0
	107	pdudt(i,k)=0.0
	108	pdvdt(i,k)=0.0
	109	pdtdt(i,k)=0.0
	110	ENDDO
	111	ENDDO
	112	c
	113	c PREPARE INPUT VARIABLES FOR ORODRAG (i.e., ORDERED FROM TOP TO BOTTOM)
	114	C CALCULATE LAYERS HEIGHT ABOVE GROUND)
	115	c
	116	DO k = 1, klev
	117	DO i = 1, klon
	118	pt(i,k) = t(i,klev-k+1)
	119	pu(i,k) = u(i,klev-k+1)
	120	pv(i,k) = v(i,klev-k+1)
	121	papmf(i,k) = pplay(i,klev-k+1)
	122	ENDDO
	123	ENDDO
	124	DO k = 1, klev+1
	125	DO i = 1, klon
	126	papmh(i,k) = paprs(i,klev-k+2)
	127	ENDDO
	128	ENDDO
	129	DO i = 1, klon
	130	zgeom(i,klev) = RD * pt(i,klev)
	131	. * LOG(papmh(i,klev+1)/papmf(i,klev))
	132	ENDDO
	133	DO k = klev-1, 1, -1
	134	DO i = 1, klon
	135	zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
	136	. * LOG(papmf(i,k+1)/papmf(i,k))
	137	ENDDO
	138	ENDDO
	139	c
	140	c CALL SSO DRAG ROUTINES
	141	c
	142	CALL orodrag_strato(klon,klev,kgwd,kdx,ktest,
	143	. dtime,
	144	. papmh, papmf, zgeom,
	145	. pt, pu, pv,
	146	. pmea, pstd, psig, pgam, pthe, ppic,pval,
	147	. pulow,pvlow,
	148	. pdudt,pdvdt,pdtdt)
	149	C
	150	C COMPUTE INCREMENTS AND STRESS FROM TENDENCIES
	151
	152	DO k = 1, klev
	153	DO i = 1, klon
	154	d_u(i,klev+1-k) = dtime*pdudt(i,k)
	155	d_v(i,klev+1-k) = dtime*pdvdt(i,k)
	156	d_t(i,klev+1-k) = dtime*pdtdt(i,k)
	157	pustr(i) = pustr(i)
	158	. +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
	159	pvstr(i) = pvstr(i)
	160	. +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
	161	ENDDO
	162	ENDDO
	163	c
	164	RETURN
	165	END
	166
	167	SUBROUTINE orodrag_strato( nlon,nlev
	168	i , kgwd, kdx, ktest
	169	r , ptsphy
	170	r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
	171	r , pmea, pstd, psig, pgam, pthe, ppic, pval
	172	c outputs
	173	r , pulow,pvlow
	174	r , pvom,pvol,pte )
	175
	176	USE dimphy
	177	IMPLICIT NONE
	178	c
	179	c
	180	c**** orodrag - does the SSO drag parametrization.
	181	c
	182	c purpose.
	183	c --------
	184	c
	185	c this routine computes the physical tendencies of the
	186	c prognostic variables u,v and t due to vertical transports by
	187	c subgridscale orographically excited gravity waves, and to
	188	c low level blocked flow drag.
	189	c
	190	c** interface.
	191	c ----------
	192	c called from drag_noro.
	193	c
	194	c the routine takes its input from the long-term storage:
	195	c u,v,t and p at t-1.
	196	c
	197	c explicit arguments :
	198	c --------------------
	199	c ==== inputs ===
	200	c nlon----input-I-Total number of horizontal points that get into physics
	201	c nlev----input-I-Number of vertical levels
	202	c
	203	c kgwd- -input-I: Total nb of points where the orography schemes are active
	204	c ktest--input-I: Flags to indicate active points
	205	c kdx----input-I: Locate the physical location of an active point.
	206	c ptsphy--input-R-Time-step (s)
	207	c paphm1--input-R: pressure at model 1/2 layer
	208	c papm1---input-R: pressure at model layer
	209	c pgeom1--input-R: Altitude of layer above ground
	210	c ptm1, pum1, pvm1--R-: t, u and v
	211	c pmea----input-R-Mean Orography (m)
	212	C pstd----input-R-SSO standard deviation (m)
	213	c psig----input-R-SSO slope
	214	c pgam----input-R-SSO Anisotropy
	215	c pthe----input-R-SSO Angle
	216	c ppic----input-R-SSO Peacks elevation (m)
	217	c pval----input-R-SSO Valleys elevation (m)
	218
	219	integer nlon,nlev,kgwd
	220	real ptsphy
	221
	222	c ==== outputs ===
	223	c pulow, pvlow -output-R: Low-level wind
	224	c
	225	c pte -----output-R: T tendency
	226	c pvom-----output-R: U tendency
	227	c pvol-----output-R: V tendency
	228	c
	229	c
	230	c Implicit Arguments:
	231	c ===================
	232	c
	233	c klon-common-I: Number of points seen by the physics
	234	c klev-common-I: Number of vertical layers
	235	c
	236	c method.
	237	c -------
	238	c
	239	c externals.
	240	c ----------
	241	integer ismin, ismax
	242	external ismin, ismax
	243	c
	244	c reference.
	245	c ----------
	246	c
	247	c author.
	248	c -------
	249	c m.miller + b.ritter e.c.m.w.f. 15/06/86.
	250	c
	251	c f.lott + m. miller e.c.m.w.f. 22/11/94
	252	c-----------------------------------------------------------------------
	253	c
	254	c
	255	cym#include "dimensions.h"
	256	cym#include "dimphy.h"
	257	#include "YOMCST.h"
	258	#include "YOEGWD.h"
	259	c-----------------------------------------------------------------------
	260	c
	261	c* 0.1 arguments
	262	c ---------
	263	c
	264	c
	265	real pte(nlon,nlev),
	266	* pvol(nlon,nlev),
	267	* pvom(nlon,nlev),
	268	* pulow(nlon),
	269	* pvlow(nlon)
	270	real pum1(nlon,nlev),
	271	* pvm1(nlon,nlev),
	272	* ptm1(nlon,nlev),
	273	* pmea(nlon),pstd(nlon),psig(nlon),
	274	* pgam(nlon),pthe(nlon),ppic(nlon),pval(nlon),
	275	* pgeom1(nlon,nlev),
	276	* papm1(nlon,nlev),
	277	* paphm1(nlon,nlev+1)
	278	c
	279	integer kdx(nlon),ktest(nlon)
	280	c-----------------------------------------------------------------------
	281	c
	282	c* 0.2 local arrays
	283	c ------------
	284	integer isect(klon),
	285	* icrit(klon),
	286	* ikcrith(klon),
	287	* ikenvh(klon),
	288	* iknu(klon),
	289	* iknu2(klon),
	290	* ikcrit(klon),
	291	* ikhlim(klon)
	292	c
	293	real ztau(klon,klev+1),
	294	* zstab(klon,klev+1),
	295	* zvph(klon,klev+1),
	296	* zrho(klon,klev+1),
	297	* zri(klon,klev+1),
	298	* zpsi(klon,klev+1),
	299	* zzdep(klon,klev)
	300	real zdudt(klon),
	301	* zdvdt(klon),
	302	* zdtdt(klon),
	303	* zdedt(klon),
	304	* zvidis(klon),
	305	* ztfr(klon),
	306	* znu(klon),
	307	* zd1(klon),
	308	* zd2(klon),
	309	* zdmod(klon)
	310
	311
	312	c local quantities:
	313
	314	integer jl,jk,ji
	315	real ztmst,zdelp,ztemp,zforc,ztend,rover
	316	real zb,zc,zconb,zabsv,zzd1,ratio,zbet,zust,zvst,zdis
	317
	318	c
	319	c------------------------------------------------------------------
	320	c
	321	c* 1. initialization
	322	c --------------
	323	c
	324	c print *,' in orodrag'
	325	100 continue
	326	c
	327	c ------------------------------------------------------------------
	328	c
	329	c* 1.1 computational constants
	330	c -----------------------
	331	c
	332	110 continue
	333	c
	334	c ztmst=twodt
	335	c if(nstep.eq.nstart) ztmst=0.5*twodt
	336	ztmst=ptsphy
	337	c ------------------------------------------------------------------
	338	c
	339	120 continue
	340	c
	341	c ------------------------------------------------------------------
	342	c
	343	c* 1.3 check whether row contains point for printing
	344	c ---------------------------------------------
	345	c
	346	130 continue
	347	c
	348	c ------------------------------------------------------------------
	349	c
	350	c* 2. precompute basic state variables.
	351	c* ---------- ----- ----- ----------
	352	c* define low level wind, project winds in plane of
	353	c* low level wind, determine sector in which to take
	354	c* the variance and set indicator for critical levels.
	355	c
	356
	357	200 continue
	358	c
	359	c
	360	c
	361	call orosetup_strato
	362	* ( nlon, nlev , ktest
	363	* , ikcrit, ikcrith, icrit, isect, ikhlim, ikenvh,iknu,iknu2
	364	* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
	365	* , zrho , zri , zstab , ztau , zvph , zpsi, zzdep
	366	* , pulow, pvlow
	367	* , pthe,pgam,pmea,ppic,pval,znu ,zd1, zd2, zdmod )
	368
	369
	370	c
	371	c
	372	c
	373	c***********************************************************
	374	c
	375	c
	376	c* 3. compute low level stresses using subcritical and
	377	c* supercritical forms.computes anisotropy coefficient
	378	c* as measure of orographic twodimensionality.
	379	c
	380	300 continue
	381	c
	382	call gwstress_strato
	383	* ( nlon , nlev
	384	* , ikcrit, isect, ikhlim, ktest, ikcrith, icrit, ikenvh, iknu
	385	* , zrho , zstab, zvph , pstd, psig, pmea, ppic, pval
	386	* , ztfr , ztau
	387	* , pgeom1,pgam,zd1,zd2,zdmod,znu)
	388
	389	c
	390	c
	391	c* 4. compute stress profile including
	392	c trapped waves, wave breaking,
	393	c linear decay in stratosphere.
	394	c
	395	400 continue
	396	c
	397	c
	398
	399	call gwprofil_strato
	400	* ( nlon , nlev
	401	* , kgwd , kdx , ktest
	402	* , ikcrit, ikcrith, icrit , ikenvh, iknu
	403	* ,iknu2 , paphm1, zrho , zstab , ztfr , zvph
	404	* , zri , ztau
	405
	406	* , zdmod , znu , psig , pgam , pstd , ppic , pval)
	407
	408	c
	409	c* 5. Compute tendencies from waves stress profile.
	410	c Compute low level blocked flow drag.
	411	c* --------------------------------------------
	412	c
	413	500 continue
	414
	415
	416	c
	417	c explicit solution at all levels for the gravity wave
	418	c implicit solution for the blocked levels
	419
	420	do 510 jl=kidia,kfdia
	421	zvidis(jl)=0.0
	422	zdudt(jl)=0.0
	423	zdvdt(jl)=0.0
	424	zdtdt(jl)=0.0
	425	510 continue
	426	c
	427
	428	do 524 jk=1,klev
	429	c
	430
	431	C WAVE STRESS
	432	C-------------
	433	c
	434	c
	435	do 523 ji=kidia,kfdia
	436
	437	if(ktest(ji).eq.1) then
	438
	439	zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
	440	ztemp=-rg(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,klev+1)zdelp)
	441
	442	zdudt(ji)=(pulow(ji)zd1(ji)-pvlow(ji)zd2(ji))*ztemp/zdmod(ji)
	443	zdvdt(ji)=(pvlow(ji)zd1(ji)+pulow(ji)zd2(ji))*ztemp/zdmod(ji)
	444	c
	445	c Control Overshoots
	446	c
	447
	448	if(jk.ge.nstra)then
	449	rover=0.10
	450	if(abs(zdudt(ji)).gt.rover*abs(pum1(ji,jk))/ztmst)
	451	C zdudt(ji)=roverabs(pum1(ji,jk))/ztmst
	452	C zdudt(ji)/(abs(zdudt(ji))+1.E-10)
	453	if(abs(zdvdt(ji)).gt.rover*abs(pvm1(ji,jk))/ztmst)
	454	C zdvdt(ji)=roverabs(pvm1(ji,jk))/ztmst
	455	C zdvdt(ji)/(abs(zdvdt(ji))+1.E-10)
	456	endif
	457
	458	rover=0.25
	459	zforc=sqrt(zdudt(ji)2+zdvdt(ji)2)
	460	ztend=sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/ztmst
	461
	462	if(zforc.ge.rover*ztend)then
	463	zdudt(ji)=roverztend/zforczdudt(ji)
	464	zdvdt(ji)=roverztend/zforczdvdt(ji)
	465	endif
	466	c
	467	c BLOCKED FLOW DRAG:
	468	C -----------------
	469	c
	470	if(jk.gt.ikenvh(ji)) then
	471	zb=1.0-0.18pgam(ji)-0.04pgam(ji)**2
	472	zc=0.48pgam(ji)+0.3pgam(ji)**2
	473	zconb=2.ztmstgkwakepsig(ji)/(4.pstd(ji))
	474	zabsv=sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/2.
	475	zzd1=zbcos(zpsi(ji,jk))2+zcsin(zpsi(ji,jk))**2
	476	ratio=(cos(zpsi(ji,jk))*2+pgam(ji)sin(zpsi(ji,jk))**2)/
	477	* (pgam(ji)cos(zpsi(ji,jk))2+sin(zpsi(ji,jk))*2)
	478	zbet=max(0.,2.-1./ratio)zconbzzdep(ji,jk)zzd1zabsv
	479	c
	480	c OPPOSED TO THE WIND
	481	c
	482	zdudt(ji)=-pum1(ji,jk)/ztmst
	483	zdvdt(ji)=-pvm1(ji,jk)/ztmst
	484	c
	485	c PERPENDICULAR TO THE SSO MAIN AXIS:
	486	C
	487	cmod zdudt(ji)=-(pum1(ji,jk)cos(pthe(ji)rpi/180.)
	488	cmod * +pvm1(ji,jk)sin(pthe(ji)rpi/180.))
	489	cmod * cos(pthe(ji)rpi/180.)/ztmst
	490	cmod zdvdt(ji)=-(pum1(ji,jk)cos(pthe(ji)rpi/180.)
	491	cmod * +pvm1(ji,jk)sin(pthe(ji)rpi/180.))
	492	cmod * sin(pthe(ji)rpi/180.)/ztmst
	493	C
	494	zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
	495	zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
	496	end if
	497	pvom(ji,jk)=zdudt(ji)
	498	pvol(ji,jk)=zdvdt(ji)
	499	zust=pum1(ji,jk)+ztmst*zdudt(ji)
	500	zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
	501	zdis=0.5(pum1(ji,jk)2+pvm1(ji,jk)2-zust2-zvst*2)
	502	zdedt(ji)=zdis/ztmst
	503	zvidis(ji)=zvidis(ji)+zdis*zdelp
	504	zdtdt(ji)=zdedt(ji)/rcpd
	505	c
	506	c NO TENDENCIES ON TEMPERATURE .....
	507	c
	508	c Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation
	509	c
	510	pte(ji,jk)=0.0
	511
	512	endif
	513
	514	523 continue
	515	524 continue
	516	c
	517	c
	518	501 continue
	519
	520	return
	521	end
	522	SUBROUTINE orosetup_strato
	523	* ( nlon , nlev , ktest
	524	* , kkcrit, kkcrith, kcrit, ksect , kkhlim
	525	* , kkenvh, kknu , kknu2
	526	* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
	527	* , prho , pri , pstab , ptau , pvph ,ppsi, pzdep
	528	* , pulow , pvlow
	529	* , ptheta, pgam, pmea, ppic, pval
	530	* , pnu , pd1 , pd2 ,pdmod )
	531	C
	532	c**** gwsetup
	533	c
	534	c purpose.
	535	c --------
	536	c SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME:
	537	C DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND
	538	C STRATIFICATION.....
	539	c
	540	c** interface.
	541	c ----------
	542	c from orodrag
	543	c
	544	c explicit arguments :
	545	c --------------------
	546	c ==== inputs ===
	547	c
	548	c nlon----input-I-Total number of horizontal points that get into physics
	549	c nlev----input-I-Number of vertical levels
	550	c ktest--input-I: Flags to indicate active points
	551	c
	552	c ptsphy--input-R-Time-step (s)
	553	c paphm1--input-R: pressure at model 1/2 layer
	554	c papm1---input-R: pressure at model layer
	555	c pgeom1--input-R: Altitude of layer above ground
	556	c ptm1, pum1, pvm1--R-: t, u and v
	557	c pmea----input-R-Mean Orography (m)
	558	C pstd----input-R-SSO standard deviation (m)
	559	c psig----input-R-SSO slope
	560	c pgam----input-R-SSO Anisotropy
	561	c pthe----input-R-SSO Angle
	562	c ppic----input-R-SSO Peacks elevation (m)
	563	c pval----input-R-SSO Valleys elevation (m)
	564
	565	c ==== outputs ===
	566	c pulow, pvlow -output-R: Low-level wind
	567	c kkcrit----I-: Security value for top of low level flow
	568	c kcrit-----I-: Critical level
	569	c ksect-----I-: Not used
	570	c kkhlim----I-: Not used
	571	c kkenvh----I-: Top of blocked flow layer
	572	c kknu------I-: Layer that sees mountain peacks
	573	c kknu2-----I-: Layer that sees mountain peacks above mountain mean
	574	c kknub-----I-: Layer that sees mountain mean above valleys
	575	c prho------R-: Density at 1/2 layers
	576	c pri-------R-: Background Richardson Number, Wind shear measured along GW stress
	577	c pstab-----R-: Brunt-Vaisala freq. at 1/2 layers
	578	c pvph------R-: Wind in plan of GW stress, Half levels.
	579	c ppsi------R-: Angle between low level wind and SS0 main axis.
	580	c pd1-------R-\| Compared the ratio of the stress
	581	c pd2-------R-\| that is along the wind to that Normal to it.
	582	c pdi define the plane of low level stress
	583	c compared to the low level wind.
	584	c see p. 108 Lott & Miller (1997).
	585	c pdmod-----R-: Norme of pdi
	586
	587	c === local arrays ===
	588	c
	589	c zvpf------R-: Wind projected in the plan of the low-level stress.
	590
	591	c ==== outputs ===
	592	c
	593	c implicit arguments : none
	594	c --------------------
	595	c
	596	c method.
	597	c -------
	598	c
	599	c
	600	c externals.
	601	c ----------
	602	c
	603	c
	604	c reference.
	605	c ----------
	606	c
	607	c see ecmwf research department documentation of the "i.f.s."
	608	c
	609	c author.
	610	c -------
	611	c
	612	c modifications.
	613	c --------------
	614	c f.lott for the new-gwdrag scheme november 1993
	615	c
	616	c-----------------------------------------------------------------------
	617	USE dimphy
	618	implicit none
	619	c
	620
	621	cym#include "dimensions.h"
	622	cym#include "dimphy.h"
	623	#include "YOMCST.h"
	624	#include "YOEGWD.h"
	625
	626	c-----------------------------------------------------------------------
	627	c
	628	c* 0.1 arguments
	629	c ---------
	630	c
	631	integer nlon,nlev
	632	integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
	633	* kkhlim(nlon),ktest(nlon),kkenvh(nlon)
	634
	635	c
	636	real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
	637	* pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
	638	* prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
	639	* ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
	640	* pzdep(nlon,klev)
	641	real pulow(nlon),pvlow(nlon),ptheta(nlon),pgam(nlon),pnu(nlon),
	642	* pd1(nlon),pd2(nlon),pdmod(nlon)
	643	real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
	644	c
	645	c-----------------------------------------------------------------------
	646	c
	647	c* 0.2 local arrays
	648	c ------------
	649	c
	650	c
	651	integer ilevh ,jl,jk
	652	real zcons1,zcons2,zhgeo,zu,zphi
	653	real zvt1,zvt2,zdwind,zwind,zdelp
	654	real zstabm,zstabp,zrhom,zrhop
	655	logical lo
	656	logical ll1(klon,klev+1)
	657	integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
	658	* kentp(klon),ncount(klon)
	659	c
	660	real zhcrit(klon,klev),zvpf(klon,klev),
	661	* zdp(klon,klev)
	662	real znorm(klon),zb(klon),zc(klon),
	663	* zulow(klon),zvlow(klon),znup(klon),znum(klon)
	664	c
	665	c ------------------------------------------------------------------
	666	c
	667	c* 1. initialization
	668	c --------------
	669	c
	670	c PRINT *,' in orosetup'
	671	100 continue
	672	c
	673	c ------------------------------------------------------------------
	674	c
	675	c* 1.1 computational constants
	676	c -----------------------
	677	c
	678	110 continue
	679	c
	680	ilevh =klev/3
	681	c
	682	zcons1=1./rd
	683	zcons2=rg**2/rcpd
	684	c
	685	c
	686	c ------------------------------------------------------------------
	687	c
	688	c* 2.
	689	c --------------
	690	c
	691	200 continue
	692	c
	693	c ------------------------------------------------------------------
	694	c
	695	c* 2.1 define low level wind, project winds in plane of
	696	c* low level wind, determine sector in which to take
	697	c* the variance and set indicator for critical levels.
	698	c
	699	c
	700	c
	701	do 2001 jl=kidia,kfdia
	702	kknu(jl) =klev
	703	kknu2(jl) =klev
	704	kknub(jl) =klev
	705	kknul(jl) =klev
	706	pgam(jl) =max(pgam(jl),gtsec)
	707	ll1(jl,klev+1)=.false.
	708	2001 continue
	709	c
	710	c Ajouter une initialisation (L. Li, le 23fev99):
	711	c
	712	do jk=klev,ilevh,-1
	713	do jl=kidia,kfdia
	714	ll1(jl,jk)= .false.
	715	ENDDO
	716	ENDDO
	717	c
	718	c* define top of low level flow
	719	c ----------------------------
	720	do 2002 jk=klev,ilevh,-1
	721	do 2003 jl=kidia,kfdia
	722	if(ktest(jl).eq.1) then
	723	lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
	724	if(lo) then
	725	kkcrit(jl)=jk
	726	endif
	727	zhcrit(jl,jk)=ppic(jl)-pval(jl)
	728	zhgeo=pgeom1(jl,jk)/rg
	729	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
	730	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
	731	kknu(jl)=jk
	732	endif
	733	if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
	734	endif
	735	2003 continue
	736	2002 continue
	737	do 2004 jk=klev,ilevh,-1
	738	do 2005 jl=kidia,kfdia
	739	if(ktest(jl).eq.1) then
	740	zhcrit(jl,jk)=ppic(jl)-pmea(jl)
	741	zhgeo=pgeom1(jl,jk)/rg
	742	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
	743	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
	744	kknu2(jl)=jk
	745	endif
	746	if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
	747	endif
	748	2005 continue
	749	2004 continue
	750	do 2006 jk=klev,ilevh,-1
	751	do 2007 jl=kidia,kfdia
	752	if(ktest(jl).eq.1) then
	753	zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
	754	zhgeo=pgeom1(jl,jk)/rg
	755	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
	756	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
	757	kknub(jl)=jk
	758	endif
	759	if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
	760	endif
	761	2007 continue
	762	2006 continue
	763	c
	764	do 2010 jl=kidia,kfdia
	765	if(ktest(jl).eq.1) then
	766	kknu(jl)=min(kknu(jl),nktopg)
	767	kknu2(jl)=min(kknu2(jl),nktopg)
	768	kknub(jl)=min(kknub(jl),nktopg)
	769	kknul(jl)=klev
	770	endif
	771	2010 continue
	772	c
	773	210 continue
	774	c
	775	cc* initialize various arrays
	776	c
	777	do 2107 jl=kidia,kfdia
	778	prho(jl,klev+1) =0.0
	779	cym correction en attendant mieux
	780	prho(jl,1) =0.0
	781	pstab(jl,klev+1) =0.0
	782	pstab(jl,1) =0.0
	783	pri(jl,klev+1) =9999.0
	784	ppsi(jl,klev+1) =0.0
	785	pri(jl,1) =0.0
	786	pvph(jl,1) =0.0
	787	pvph(jl,klev+1) =0.0
	788	cym correction en attendant mieux
	789	cym pvph(jl,klev) =0.0
	790	pulow(jl) =0.0
	791	pvlow(jl) =0.0
	792	zulow(jl) =0.0
	793	zvlow(jl) =0.0
	794	kkcrith(jl) =klev
	795	kkenvh(jl) =klev
	796	kentp(jl) =klev
	797	kcrit(jl) =1
	798	ncount(jl) =0
	799	ll1(jl,klev+1) =.false.
	800	2107 continue
	801	c
	802	c* define flow density and stratification (rho and N2)
	803	c at semi layers.
	804	c -------------------------------------------------------
	805	c
	806	do 223 jk=klev,2,-1
	807	do 222 jl=kidia,kfdia
	808	if(ktest(jl).eq.1) then
	809	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
	810	prho(jl,jk)=2.paphm1(jl,jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
	811	pstab(jl,jk)=2.zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))
	812	* (1.-rcpdprho(jl,jk)(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
	813	pstab(jl,jk)=max(pstab(jl,jk),gssec)
	814	endif
	815	222 continue
	816	223 continue
	817	c
	818	c********************************************************************
	819	c
	820	c* define Low level flow (between ground and peacks-valleys)
	821	c ---------------------------------------------------------
	822	do 2115 jk=klev,ilevh,-1
	823	do 2116 jl=kidia,kfdia
	824	if(ktest(jl).eq.1) then
	825	if(jk.ge.kknu2(jl).and.jk.le.kknul(jl)) then
	826	pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	827	pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	828	pstab(jl,klev+1)=pstab(jl,klev+1)
	829	c +pstab(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	830	prho(jl,klev+1)=prho(jl,klev+1)
	831	c +prho(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	832	end if
	833	endif
	834	2116 continue
	835	2115 continue
	836	do 2110 jl=kidia,kfdia
	837	if(ktest(jl).eq.1) then
	838	pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
	839	pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
	840	znorm(jl)=max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
	841	pvph(jl,klev+1)=znorm(jl)
	842	pstab(jl,klev+1)=pstab(jl,klev+1)
	843	c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
	844	prho(jl,klev+1)=prho(jl,klev+1)
	845	c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
	846	endif
	847	2110 continue
	848
	849	c
	850	c******* setup orography orientation relative to the low level
	851	C wind and define parameters of the Anisotropic wave stress.
	852	c
	853	do 2112 jl=kidia,kfdia
	854	if(ktest(jl).eq.1) then
	855	lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
	856	if(lo) then
	857	zu=pulow(jl)+2.*gvsec
	858	else
	859	zu=pulow(jl)
	860	endif
	861	zphi=atan(pvlow(jl)/zu)
	862	ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
	863	zb(jl)=1.-0.18pgam(jl)-0.04pgam(jl)**2
	864	zc(jl)=0.48pgam(jl)+0.3pgam(jl)**2
	865	pd1(jl)=zb(jl)-(zb(jl)-zc(jl))(sin(ppsi(jl,klev+1))*2)
	866	pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))
	867	* *cos(ppsi(jl,klev+1))
	868	pdmod(jl)=sqrt(pd1(jl)2+pd2(jl)2)
	869	endif
	870	2112 continue
	871	c
	872	c ********** projet flow in plane of lowlevel stress ***********
	873	C ********** Find critical levels... ***********
	874	c
	875	do 213 jk=1,klev
	876	do 212 jl=kidia,kfdia
	877	if(ktest(jl).eq.1) then
	878	zvt1 =pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk)
	879	zvt2 =-pvlow(jl)pum1(jl,jk)+pulow(jl)pvm1(jl,jk)
	880	zvpf(jl,jk)=(zvt1pd1(jl)+zvt2pd2(jl))/(znorm(jl)*pdmod(jl))
	881	endif
	882	ptau(jl,jk) =0.0
	883	pzdep(jl,jk) =0.0
	884	ppsi(jl,jk) =0.0
	885	ll1(jl,jk) =.false.
	886	212 continue
	887	213 continue
	888	do 215 jk=2,klev
	889	do 214 jl=kidia,kfdia
	890	if(ktest(jl).eq.1) then
	891	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
	892	pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
	893	* (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
	894	* /zdp(jl,jk)
	895	if(pvph(jl,jk).lt.gvsec) then
	896	pvph(jl,jk)=gvsec
	897	kcrit(jl)=jk
	898	endif
	899	endif
	900	214 continue
	901	215 continue
	902	c
	903	c* 2.3 mean flow richardson number.
	904	c
	905	230 continue
	906	c
	907	do 232 jk=2,klev
	908	do 231 jl=kidia,kfdia
	909	if(ktest(jl).eq.1) then
	910	zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
	911	pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
	912	* /(rgprho(jl,jk)zdwind))**2
	913	pri(jl,jk)=max(pri(jl,jk),grcrit)
	914	endif
	915	231 continue
	916	232 continue
	917
	918	c
	919	c
	920	c* define top of 'envelope' layer
	921	c ----------------------------
	922
	923	do 233 jl=kidia,kfdia
	924	pnu (jl)=0.0
	925	znum(jl)=0.0
	926	233 continue
	927
	928	do 234 jk=2,klev-1
	929	do 234 jl=kidia,kfdia
	930
	931	if(ktest(jl).eq.1) then
	932
	933	if (jk.ge.kknu2(jl)) then
	934
	935	znum(jl)=pnu(jl)
	936	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
	937	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
	938	zwind=max(sqrt(zwind**2),gvsec)
	939	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
	940	zstabm=sqrt(max(pstab(jl,jk ),gssec))
	941	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
	942	zrhom=prho(jl,jk )
	943	zrhop=prho(jl,jk+1)
	944	pnu(jl) = pnu(jl) + (zdelp/rg)*
	945	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
	946	if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
	947	* .and.(kkenvh(jl).eq.klev))
	948	* kkenvh(jl)=jk
	949
	950	endif
	951
	952	endif
	953
	954	234 continue
	955
	956	c calculation of a dynamical mixing height for when the waves
	957	C BREAK AT LOW LEVEL: The drag will be repartited over
	958	C a depths that depends on waves vertical wavelength,
	959	C not just between two adjacent model layers.
	960	c of gravity waves:
	961
	962	do 235 jl=kidia,kfdia
	963	znup(jl)=0.0
	964	znum(jl)=0.0
	965	235 continue
	966
	967	do 236 jk=klev-1,2,-1
	968	do 236 jl=kidia,kfdia
	969
	970	if(ktest(jl).eq.1) then
	971
	972	znum(jl)=znup(jl)
	973	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
	974	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
	975	zwind=max(sqrt(zwind**2),gvsec)
	976	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
	977	zstabm=sqrt(max(pstab(jl,jk ),gssec))
	978	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
	979	zrhom=prho(jl,jk )
	980	zrhop=prho(jl,jk+1)
	981	znup(jl) = znup(jl) + (zdelp/rg)*
	982	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
	983	if((znum(jl).le.rpi/4.).and.(znup(jl).gt.rpi/4.)
	984	* .and.(kkcrith(jl).eq.klev))
	985	* kkcrith(jl)=jk
	986
	987	endif
	988
	989	236 continue
	990
	991	do 237 jl=kidia,kfdia
	992	if(ktest(jl).eq.1) then
	993	kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
	994	kkcrith(jl)=max0(kkcrith(jl),kknu(jl))
	995	if(kcrit(jl).ge.kkcrith(jl))kcrit(jl)=1
	996	endif
	997	237 continue
	998	c
	999	c directional info for flow blocking *************************
	1000	c
	1001	do 251 jk=1,klev
	1002	do 252 jl=kidia,kfdia
	1003	if(ktest(jl).eq.1) then
	1004	lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
	1005	if(lo) then
	1006	zu=pum1(jl,jk)+2.*gvsec
	1007	else
	1008	zu=pum1(jl,jk)
	1009	endif
	1010	zphi=atan(pvm1(jl,jk)/zu)
	1011	ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
	1012	endif
	1013	252 continue
	1014	251 continue
	1015
	1016	c forms the vertical 'leakiness' **************************
	1017
	1018	do 254 jk=ilevh,klev
	1019	do 253 jl=kidia,kfdia
	1020	if(ktest(jl).eq.1) then
	1021	pzdep(jl,jk)=0
	1022	if(jk.ge.kkenvh(jl).and.kkenvh(jl).ne.klev) then
	1023	pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl) )-pgeom1(jl, jk))/
	1024	* (pgeom1(jl,kkenvh(jl) )-pgeom1(jl,klev))
	1025	end if
	1026	endif
	1027	253 continue
	1028	254 continue
	1029
	1030	return
	1031	end
	1032	SUBROUTINE gwstress_strato
	1033	* ( nlon , nlev
	1034	* , kkcrit, ksect, kkhlim, ktest, kkcrith, kcrit, kkenvh
	1035	* , kknu
	1036	* , prho , pstab , pvph , pstd, psig
	1037	* , pmea , ppic , pval , ptfr , ptau
	1038	* , pgeom1 , pgamma , pd1 , pd2 , pdmod , pnu )
	1039	c
	1040	c**** gwstress
	1041	c
	1042	c purpose.
	1043	c --------
	1044	c Compute the surface stress due to Gravity Waves, according
	1045	c to the Phillips (1979) theory of 3-D flow above
	1046	c anisotropic elliptic ridges.
	1047
	1048	C The stress is reduced two account for cut-off flow over
	1049	C hill. The flow only see that part of the ridge located
	1050	c above the blocked layer (see zeff).
	1051	c
	1052	c** interface.
	1053	c ----------
	1054	c call gwstress from gwdrag
	1055	c
	1056	c explicit arguments :
	1057	c --------------------
	1058	c ==== inputs ===
	1059	c ==== outputs ===
	1060	c
	1061	c implicit arguments : none
	1062	c --------------------
	1063	c
	1064	c method.
	1065	c -------
	1066	c
	1067	c
	1068	c externals.
	1069	c ----------
	1070	c
	1071	c
	1072	c reference.
	1073	c ----------
	1074	c
	1075	c LOTT and MILLER (1997) & LOTT (1999)
	1076	c
	1077	c author.
	1078	c -------
	1079	c
	1080	c modifications.
	1081	c --------------
	1082	c f. lott put the new gwd on ifs 22/11/93
	1083	c
	1084	c-----------------------------------------------------------------------
	1085	USE dimphy
	1086	implicit none
	1087
	1088	cym#include "dimensions.h"
	1089	cym#include "dimphy.h"
	1090	#include "YOMCST.h"
	1091	#include "YOEGWD.h"
	1092
	1093	c-----------------------------------------------------------------------
	1094	c
	1095	c* 0.1 arguments
	1096	c ---------
	1097	c
	1098	integer nlon,nlev
	1099	integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
	1100	* kkhlim(nlon),ktest(nlon),kkenvh(nlon),kknu(nlon)
	1101	c
	1102	real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1),
	1103	* pvph(nlon,nlev+1),ptfr(nlon),
	1104	* pgeom1(nlon,nlev),pstd(nlon)
	1105	c
	1106	real pd1(nlon),pd2(nlon),pnu(nlon),psig(nlon),pgamma(nlon)
	1107	real pmea(nlon),ppic(nlon),pval(nlon)
	1108	real pdmod(nlon)
	1109	c
	1110	c-----------------------------------------------------------------------
	1111	c
	1112	c* 0.2 local arrays
	1113	c ------------
	1114	c zeff--real: effective height seen by the flow when there is blocking
	1115
	1116	integer jl
	1117	real zeff
	1118	c
	1119	c-----------------------------------------------------------------------
	1120	c
	1121	c* 0.3 functions
	1122	c ---------
	1123	c ------------------------------------------------------------------
	1124	c
	1125	c* 1. initialization
	1126	c --------------
	1127	c
	1128	c PRINT *,' in gwstress'
	1129	100 continue
	1130	c
	1131	c* 3.1 gravity wave stress.
	1132	c
	1133	300 continue
	1134	c
	1135	c
	1136	do 301 jl=kidia,kfdia
	1137	if(ktest(jl).eq.1) then
	1138
	1139	c effective mountain height above the blocked flow
	1140
	1141	zeff=ppic(jl)-pval(jl)
	1142	if(kkenvh(jl).lt.klev)then
	1143	zeff=amin1(GFRCRIT*pvph(jl,klev+1)/sqrt(pstab(jl,klev+1))
	1144	c ,zeff)
	1145	endif
	1146
	1147
	1148	ptau(jl,klev+1)=gkdrag*prho(jl,klev+1)
	1149	* psig(jl)pdmod(jl)/4./pstd(jl)
	1150	* pvph(jl,klev+1)sqrt(pstab(jl,klev+1))
	1151	* zeff*2
	1152
	1153
	1154	c too small value of stress or low level flow include critical level
	1155	c or low level flow: gravity wave stress nul.
	1156
	1157	c lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl))
	1158	c * .or.(pvph(jl,klev+1).lt.gvcrit)
	1159	c if(lo) ptau(jl,klev+1)=0.0
	1160
	1161	c print *,jl,ptau(jl,klev+1)
	1162
	1163	else
	1164
	1165	ptau(jl,klev+1)=0.0
	1166
	1167	endif
	1168
	1169	301 continue
	1170
	1171	c write(21)(ptau(jl,klev+1),jl=kidia,kfdia)
	1172
	1173	return
	1174	end
	1175
	1176	subroutine gwprofil_strato
	1177	* ( nlon, nlev
	1178	* , kgwd ,kdx , ktest
	1179	* , kkcrit, kkcrith, kcrit , kkenvh, kknu,kknu2
	1180	* , paphm1, prho , pstab , ptfr , pvph , pri , ptau
	1181	* , pdmod , pnu , psig ,pgamma, pstd, ppic,pval)
	1182
	1183	C**** gwprofil
	1184	C
	1185	C purpose.
	1186	C --------
	1187	C
	1188	C** interface.
	1189	C ----------
	1190	C from gwdrag
	1191	C
	1192	C explicit arguments :
	1193	C --------------------
	1194	C ==== inputs ===
	1195	C
	1196	C ==== outputs ===
	1197	C
	1198	C implicit arguments : none
	1199	C --------------------
	1200	C
	1201	C method:
	1202	C -------
	1203	C the stress profile for gravity waves is computed as follows:
	1204	C it decreases linearly with heights from the ground
	1205	C to the low-level indicated by kkcrith,
	1206	C to simulates lee waves or
	1207	C low-level gravity wave breaking.
	1208	C above it is constant, except when the waves encounter a critical
	1209	C level (kcrit) or when they break.
	1210	C The stress is also uniformly distributed above the level
	1211	C nstra.
	1212	C
	1213	USE dimphy
	1214	IMPLICIT NONE
	1215
	1216	cym#include "dimensions.h"
	1217	cym#include "dimphy.h"
	1218	#include "YOMCST.h"
	1219	#include "YOEGWD.h"
	1220
	1221	C-----------------------------------------------------------------------
	1222	C
	1223	C* 0.1 ARGUMENTS
	1224	C ---------
	1225	C
	1226	integer nlon,nlev,kgwd
	1227	integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon)
	1228	* ,kdx(nlon),ktest(nlon)
	1229	* ,kkenvh(nlon),kknu(nlon),kknu2(nlon)
	1230	C
	1231	real paphm1(nlon,nlev+1), pstab(nlon,nlev+1),
	1232	* prho (nlon,nlev+1), pvph (nlon,nlev+1),
	1233	* pri (nlon,nlev+1), ptfr (nlon), ptau(nlon,nlev+1)
	1234
	1235	real pdmod (nlon) , pnu (nlon) , psig(nlon),
	1236	* pgamma(nlon) , pstd(nlon) , ppic(nlon), pval(nlon)
	1237
	1238	C-----------------------------------------------------------------------
	1239	C
	1240	C* 0.2 local arrays
	1241	C ------------
	1242	C
	1243	integer jl,jk
	1244	real zsqr,zalfa,zriw,zdel,zb,zalpha,zdz2n,zdelp,zdelpt
	1245
	1246	real zdz2 (klon,klev) , znorm(klon) , zoro(klon)
	1247	real ztau (klon,klev+1)
	1248	C
	1249	C-----------------------------------------------------------------------
	1250	C
	1251	C* 1. INITIALIZATION
	1252	C --------------
	1253	C
	1254	C print *,' entree gwprofil'
	1255	100 CONTINUE
	1256	C
	1257	C
	1258	C* COMPUTATIONAL CONSTANTS.
	1259	C ------------- ----------
	1260	C
	1261	do 400 jl=kidia,kfdia
	1262	if(ktest(jl).eq.1)then
	1263	zoro(jl)=psig(jl)*pdmod(jl)/4./pstd(jl)
	1264	ztau(jl,klev+1)=ptau(jl,klev+1)
	1265	c print *,jl,ptau(jl,klev+1)
	1266	ztau(jl,kkcrith(jl))=grahilo*ptau(jl,klev+1)
	1267	endif
	1268	400 continue
	1269
	1270	C
	1271	do 430 jk=klev+1,1,-1
	1272	C
	1273	C
	1274	C* 4.1 constant shear stress until top of the
	1275	C low-level breaking/trapped layer
	1276	410 CONTINUE
	1277	C
	1278	do 411 jl=kidia,kfdia
	1279	if(ktest(jl).eq.1)then
	1280	if(jk.gt.kkcrith(jl)) then
	1281	zdelp=paphm1(jl,jk)-paphm1(jl,klev+1)
	1282	zdelpt=paphm1(jl,kkcrith(jl))-paphm1(jl,klev+1)
	1283	ptau(jl,jk)=ztau(jl,klev+1)+zdelp/zdelpt*
	1284	c (ztau(jl,kkcrith(jl))-ztau(jl,klev+1))
	1285	else
	1286	ptau(jl,jk)=ztau(jl,kkcrith(jl))
	1287	endif
	1288	endif
	1289	411 continue
	1290	C
	1291	C* 4.15 constant shear stress until the top of the
	1292	C low level flow layer.
	1293	415 continue
	1294	C
	1295	C
	1296	C* 4.2 wave displacement at next level.
	1297	C
	1298	420 continue
	1299	C
	1300	430 continue
	1301
	1302	C
	1303	C* 4.4 wave richardson number, new wave displacement
	1304	C* and stress: breaking evaluation and critical
	1305	C level
	1306	C
	1307
	1308	do 440 jk=klev,1,-1
	1309
	1310	do 441 jl=kidia,kfdia
	1311	if(ktest(jl).eq.1)then
	1312	znorm(jl)=prho(jl,jk)sqrt(pstab(jl,jk))pvph(jl,jk)
	1313	zdz2(jl,jk)=ptau(jl,jk)/amax1(znorm(jl),gssec)/zoro(jl)
	1314	endif
	1315	441 continue
	1316
	1317	do 442 jl=kidia,kfdia
	1318	if(ktest(jl).eq.1)then
	1319	if(jk.lt.kkcrith(jl)) then
	1320	if((ptau(jl,jk+1).lt.gtsec).or.(jk.le.kcrit(jl))) then
	1321	ptau(jl,jk)=0.0
	1322	else
	1323	zsqr=sqrt(pri(jl,jk))
	1324	zalfa=sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl,jk)
	1325	zriw=pri(jl,jk)(1.-zalfa)/(1+zalfazsqr)**2
	1326	if(zriw.lt.grcrit) then
	1327	C print *,' breaking!!!',ptau(jl,jk)
	1328	zdel=4./zsqr/grcrit+1./grcrit**2+4./grcrit
	1329	zb=1./grcrit+2./zsqr
	1330	zalpha=0.5*(-zb+sqrt(zdel))
	1331	zdz2n=(pvph(jl,jk)zalpha)*2/pstab(jl,jk)
	1332	ptau(jl,jk)=znorm(jl)zdz2nzoro(jl)
	1333	endif
	1334
	1335	ptau(jl,jk)=amin1(ptau(jl,jk),ptau(jl,jk+1))
	1336
	1337	endif
	1338	endif
	1339	endif
	1340	442 continue
	1341	440 continue
	1342
	1343	C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
	1344
	1345	do 530 jl=kidia,kfdia
	1346	if(ktest(jl).eq.1)then
	1347	ztau(jl,kkcrith(jl)-1)=ptau(jl,kkcrith(jl)-1)
	1348	ztau(jl,nstra)=ptau(jl,nstra)
	1349	endif
	1350	530 continue
	1351
	1352	do 531 jk=1,klev
	1353
	1354	do 532 jl=kidia,kfdia
	1355	if(ktest(jl).eq.1)then
	1356
	1357	if(jk.gt.kkcrith(jl)-1)then
	1358
	1359	zdelp=paphm1(jl,jk)-paphm1(jl,klev+1 )
	1360	zdelpt=paphm1(jl,kkcrith(jl)-1)-paphm1(jl,klev+1 )
	1361	ptau(jl,jk)=ztau(jl,klev+1 ) +
	1362	. (ztau(jl,kkcrith(jl)-1)-ztau(jl,klev+1 ) )*
	1363	. zdelp/zdelpt
	1364
	1365	endif
	1366	endif
	1367
	1368	532 continue
	1369
	1370	C REORGANISATION AT THE MODEL TOP....
	1371
	1372	do 533 jl=kidia,kfdia
	1373	if(ktest(jl).eq.1)then
	1374
	1375	if(jk.lt.nstra)then
	1376
	1377	zdelp =paphm1(jl,nstra)
	1378	zdelpt=paphm1(jl,jk)
	1379	ptau(jl,jk)=ztau(jl,nstra)*zdelpt/zdelp
	1380	c ptau(jl,jk)=ztau(jl,nstra)
	1381
	1382	endif
	1383
	1384	endif
	1385
	1386	533 continue
	1387
	1388
	1389	531 continue
	1390
	1391
	1392	123 format(i4,1x,20(f6.3,1x))
	1393
	1394
	1395	return
	1396	end
	1397	subroutine lift_noro_strato (nlon,nlev,dtime,paprs,pplay,
	1398	i plat,pmea,pstd, psig, pgam, pthe, ppic,pval,
	1399	i kgwd,kdx,ktest,
	1400	i t, u, v,
	1401	o pulow, pvlow, pustr, pvstr,
	1402	o d_t, d_u, d_v)
	1403	c
	1404	USE dimphy
	1405	implicit none
	1406	c======================================================================
	1407	c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
	1408	c Object: Mountain lift interface (enhanced vortex stretching).
	1409	c Made necessary because:
	1410	C 1. in the LMD-GCM Layers are from bottom to top,
	1411	C contrary to most European GCM.
	1412	c 2. the altitude above ground of each model layers
	1413	c needs to be known (variable zgeom)
	1414	c======================================================================
	1415	c Explicit Arguments:
	1416	c ==================
	1417	c nlon----input-I-Total number of horizontal points that get into physics
	1418	c nlev----input-I-Number of vertical levels
	1419	c dtime---input-R-Time-step (s)
	1420	c paprs---input-R-Pressure in semi layers (Pa)
	1421	c pplay---input-R-Pressure model-layers (Pa)
	1422	c t-------input-R-temperature (K)
	1423	c u-------input-R-Horizontal wind (m/s)
	1424	c v-------input-R-Meridional wind (m/s)
	1425	c pmea----input-R-Mean Orography (m)
	1426	C pstd----input-R-SSO standard deviation (m)
	1427	c psig----input-R-SSO slope
	1428	c pgam----input-R-SSO Anisotropy
	1429	c pthe----input-R-SSO Angle
	1430	c ppic----input-R-SSO Peacks elevation (m)
	1431	c pval----input-R-SSO Valleys elevation (m)
	1432	c
	1433	c kgwd- -input-I: Total nb of points where the orography schemes are active
	1434	c ktest--input-I: Flags to indicate active points
	1435	c kdx----input-I: Locate the physical location of an active point.
	1436
	1437	c pulow, pvlow -output-R: Low-level wind
	1438	c pustr, pvstr -output-R: Surface stress due to SSO drag (Pa)
	1439	c
	1440	c d_t-----output-R: T increment
	1441	c d_u-----output-R: U increment
	1442	c d_v-----output-R: V increment
	1443	c
	1444	c Implicit Arguments:
	1445	c ===================
	1446	c
	1447	c iim--common-I: Number of longitude intervals
	1448	c jjm--common-I: Number of latitude intervals
	1449	c klon-common-I: Number of points seen by the physics
	1450	c (iim+1)*(jjm+1) for instance
	1451	c klev-common-I: Number of vertical layers
	1452	c======================================================================
	1453	c Local Variables:
	1454	c ================
	1455	c
	1456	c zgeom-----R: Altitude of layer above ground
	1457	c pt, pu, pv --R: t u v from top to bottom
	1458	c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
	1459	c papmf: pressure at model layer (from top to bottom)
	1460	c papmh: pressure at model 1/2 layer (from top to bottom)
	1461	c
	1462	c======================================================================
	1463
	1464	cym#include "dimensions.h"
	1465	cym#include "dimphy.h"
	1466	#include "YOMCST.h"
	1467	#include "YOEGWD.h"
	1468	c
	1469	c ARGUMENTS
	1470	c
	1471	INTEGER nlon,nlev
	1472	REAL dtime
	1473	REAL paprs(klon,klev+1)
	1474	REAL pplay(klon,klev)
	1475	REAL plat(nlon),pmea(nlon)
	1476	REAL pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
	1477	REAL ppic(nlon),pval(nlon)
	1478	REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
	1479	REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
	1480	REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
	1481	c
	1482	INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
	1483	c
	1484	c Variables locales:
	1485	c
	1486	REAL zgeom(klon,klev)
	1487	REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
	1488	REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
	1489	REAL papmf(klon,klev),papmh(klon,klev+1)
	1490	c
	1491	c initialiser les variables de sortie (pour securite)
	1492	c
	1493
	1494	c print *,'in lift_noro'
	1495	DO i = 1,klon
	1496	pulow(i) = 0.0
	1497	pvlow(i) = 0.0
	1498	pustr(i) = 0.0
	1499	pvstr(i) = 0.0
	1500	ENDDO
	1501	DO k = 1, klev
	1502	DO i = 1, klon
	1503	d_t(i,k) = 0.0
	1504	d_u(i,k) = 0.0
	1505	d_v(i,k) = 0.0
	1506	pdudt(i,k)=0.0
	1507	pdvdt(i,k)=0.0
	1508	pdtdt(i,k)=0.0
	1509	ENDDO
	1510	ENDDO
	1511	c
	1512	c preparer les variables d'entree (attention: l'ordre des niveaux
	1513	c verticaux augmente du haut vers le bas)
	1514	c
	1515	DO k = 1, klev
	1516	DO i = 1, klon
	1517	pt(i,k) = t(i,klev-k+1)
	1518	pu(i,k) = u(i,klev-k+1)
	1519	pv(i,k) = v(i,klev-k+1)
	1520	papmf(i,k) = pplay(i,klev-k+1)
	1521	ENDDO
	1522	ENDDO
	1523	DO k = 1, klev+1
	1524	DO i = 1, klon
	1525	papmh(i,k) = paprs(i,klev-k+2)
	1526	ENDDO
	1527	ENDDO
	1528	DO i = 1, klon
	1529	zgeom(i,klev) = RD * pt(i,klev)
	1530	. * LOG(papmh(i,klev+1)/papmf(i,klev))
	1531	ENDDO
	1532	DO k = klev-1, 1, -1
	1533	DO i = 1, klon
	1534	zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
	1535	. * LOG(papmf(i,k+1)/papmf(i,k))
	1536	ENDDO
	1537	ENDDO
	1538	c
	1539	c appeler la routine principale
	1540	c
	1541
	1542	CALL OROLIFT_strato(klon,klev,kgwd,kdx,ktest,
	1543	. dtime,
	1544	. papmh, papmf, zgeom,
	1545	. pt, pu, pv,
	1546	. plat,pmea, pstd, psig, pgam, pthe, ppic,pval,
	1547	. pulow,pvlow,
	1548	. pdudt,pdvdt,pdtdt)
	1549	C
	1550	DO k = 1, klev
	1551	DO i = 1, klon
	1552	d_u(i,klev+1-k) = dtime*pdudt(i,k)
	1553	d_v(i,klev+1-k) = dtime*pdvdt(i,k)
	1554	d_t(i,klev+1-k) = dtime*pdtdt(i,k)
	1555	pustr(i) = pustr(i)
	1556	. +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
	1557	pvstr(i) = pvstr(i)
	1558	. +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
	1559	ENDDO
	1560	ENDDO
	1561
	1562	c print *,' out lift_noro'
	1563	c
	1564	RETURN
	1565	END
	1566	subroutine orolift_strato( nlon,nlev
	1567	I , kgwd, kdx, ktest
	1568	R , ptsphy
	1569	R , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
	1570	R , plat
	1571	R , pmea, pstd, psig, pgam, pthe,ppic,pval
	1572	C OUTPUTS
	1573	R , pulow,pvlow
	1574	R , pvom,pvol,pte )
	1575
	1576	C
	1577	C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
	1578	C
	1579	C PURPOSE.
	1580	C --------
	1581	C this routine computes the physical tendencies of the
	1582	C prognostic variables u,v when enhanced vortex stretching
	1583	C is needed.
	1584	C
	1585	C** INTERFACE.
	1586	C ----------
	1587	C CALLED FROM *lift_noro
	1588	c explicit arguments :
	1589	c --------------------
	1590	c ==== inputs ===
	1591	c nlon----input-I-Total number of horizontal points that get into physics
	1592	c nlev----input-I-Number of vertical levels
	1593	c
	1594	c kgwd- -input-I: Total nb of points where the orography schemes are active
	1595	c ktest--input-I: Flags to indicate active points
	1596	c kdx----input-I: Locate the physical location of an active point.
	1597	c ptsphy--input-R-Time-step (s)
	1598	c paphm1--input-R: pressure at model 1/2 layer
	1599	c papm1---input-R: pressure at model layer
	1600	c pgeom1--input-R: Altitude of layer above ground
	1601	c ptm1, pum1, pvm1--R-: t, u and v
	1602	c pmea----input-R-Mean Orography (m)
	1603	C pstd----input-R-SSO standard deviation (m)
	1604	c psig----input-R-SSO slope
	1605	c pgam----input-R-SSO Anisotropy
	1606	c pthe----input-R-SSO Angle
	1607	c ppic----input-R-SSO Peacks elevation (m)
	1608	c pval----input-R-SSO Valleys elevation (m)
	1609	c plat----input-R-Latitude (degree)
	1610	c
	1611	c ==== outputs ===
	1612	c pulow, pvlow -output-R: Low-level wind
	1613	c
	1614	c pte -----output-R: T tendency
	1615	c pvom-----output-R: U tendency
	1616	c pvol-----output-R: V tendency
	1617	c
	1618	c
	1619	c Implicit Arguments:
	1620	c ===================
	1621	c
	1622	c klon-common-I: Number of points seen by the physics
	1623	c klev-common-I: Number of vertical layers
	1624	c
	1625
	1626	C ----------
	1627	C
	1628	C AUTHOR.
	1629	C -------
	1630	C F.LOTT LMD 22/11/95
	1631	C
	1632	USE dimphy
	1633	implicit none
	1634	C
	1635	C
	1636	cym#include "dimensions.h"
	1637	cym#include "dimphy.h"
	1638	#include "YOMCST.h"
	1639	#include "YOEGWD.h"
	1640	C-----------------------------------------------------------------------
	1641	C
	1642	C* 0.1 ARGUMENTS
	1643	C ---------
	1644	C
	1645	C
	1646	integer nlon,nlev,kgwd
	1647	real ptsphy
	1648	real pte(nlon,nlev),
	1649	* pvol(nlon,nlev),
	1650	* pvom(nlon,nlev),
	1651	* pulow(nlon),
	1652	* pvlow(nlon)
	1653	real pum1(nlon,nlev),
	1654	* pvm1(nlon,nlev),
	1655	* ptm1(nlon,nlev),
	1656	* plat(nlon),pmea(nlon),
	1657	* pstd(nlon),psig(nlon),pgam(nlon),
	1658	* pthe(nlon),ppic(nlon),pval(nlon),
	1659	* pgeom1(nlon,nlev),
	1660	* papm1(nlon,nlev),
	1661	* paphm1(nlon,nlev+1)
	1662	C
	1663	INTEGER KDX(NLON),KTEST(NLON)
	1664	C-----------------------------------------------------------------------
	1665	C
	1666	C* 0.2 local arrays
	1667
	1668	integer jl,ilevh,jk
	1669	real zhgeo,zdelp,zslow,zsqua,zscav,zbet
	1670	C ------------
	1671	integer iknub(klon),
	1672	* iknul(klon)
	1673	logical ll1(klon,klev+1)
	1674	C
	1675	real ztau(klon,klev+1),
	1676	* ztav(klon,klev+1),
	1677	* zrho(klon,klev+1)
	1678	real zdudt(klon),
	1679	* zdvdt(klon)
	1680	real zhcrit(klon,klev)
	1681
	1682	logical lifthigh
	1683	real zcons1,ztmst
[1403]	1684	CHARACTER (LEN=20) :: modname='orolift_strato'
	1685	CHARACTER (LEN=80) :: abort_message
[1001]	1686
[1403]	1687
[1001]	1688	C-----------------------------------------------------------------------
	1689	C
	1690	C* 1.1 initialisations
	1691	C ---------------
	1692
	1693	lifthigh=.false.
	1694
[1403]	1695	if(nlon.ne.klon.or.nlev.ne.klev) then
	1696	abort_message = 'pb dimension'
	1697	CALL abort_gcm (modname,abort_message,1)
	1698	ENDIF
[1001]	1699	zcons1=1./rd
	1700	ztmst=ptsphy
	1701	C
	1702	do 1001 jl=kidia,kfdia
	1703	zrho(jl,klev+1) =0.0
	1704	pulow(jl) =0.0
	1705	pvlow(jl) =0.0
	1706	iknub(JL) =klev
	1707	iknul(JL) =klev
	1708	ilevh=klev/3
	1709	ll1(jl,klev+1)=.false.
	1710	do 1000 jk=1,klev
	1711	pvom(jl,jk)=0.0
	1712	pvol(jl,jk)=0.0
	1713	pte (jl,jk)=0.0
	1714	1000 continue
	1715	1001 continue
	1716
	1717	C
	1718	C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
	1719	C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
	1720	C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
	1721	C
	1722	C
	1723	C
	1724	do 2006 jk=klev,1,-1
	1725	do 2007 jl=kidia,kfdia
	1726	if(ktest(jl).eq.1) then
	1727	zhcrit(jl,jk)=amax1(ppic(jl)-pval(jl),100.)
	1728	zhgeo=pgeom1(jl,jk)/rg
	1729	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
	1730	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
	1731	iknub(jl)=jk
	1732	endif
	1733	endif
	1734	2007 continue
	1735	2006 continue
	1736	C
	1737
	1738	do 2010 jl=kidia,kfdia
	1739	if(ktest(jl).eq.1) then
	1740	iknub(jl)=max(iknub(jl),klev/2)
	1741	iknul(jl)=max(iknul(jl),2*klev/3)
	1742	if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
	1743	if(iknub(jl).eq.nktopg) iknul(jl)=klev
	1744	if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
	1745	endif
	1746	2010 continue
	1747
	1748	do 223 jk=klev,2,-1
	1749	do 222 jl=kidia,kfdia
	1750	zrho(jl,jk)=2.paphm1(jl,jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
	1751	222 continue
	1752	223 continue
	1753	c print *,' dans orolift: 223'
	1754
	1755	C********************************************************************
	1756	C
	1757	c* define low level flow
	1758	C -------------------
	1759	do 2115 jk=klev,1,-1
	1760	do 2116 jl=kidia,kfdia
	1761	if(ktest(jl).eq.1) THEN
	1762	if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
	1763	pulow(JL)=pulow(JL)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	1764	pvlow(JL)=pvlow(JL)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	1765	zrho(JL,klev+1)=zrho(JL,klev+1)
	1766	* +zrho(JL,JK)*(paphm1(jl,jk+1)-paphm1(jl,jk))
	1767	endif
	1768	endif
	1769	2116 continue
	1770	2115 continue
	1771	do 2110 jl=kidia,kfdia
	1772	if(ktest(jl).eq.1) then
	1773	pulow(JL)=pulow(JL)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
	1774	pvlow(JL)=pvlow(JL)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
	1775	zrho(JL,klev+1)=zrho(Jl,klev+1)
	1776	* /(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
	1777	endif
	1778	2110 continue
	1779
	1780
	1781	200 continue
	1782
	1783	C***********************************************************
	1784	C
	1785	C* 3. COMPUTE MOUNTAIN LIFT
	1786	C
	1787	300 continue
	1788	C
	1789	do 301 jl=kidia,kfdia
	1790	if(ktest(jl).eq.1) then
	1791	ztau(jl,klev+1)= - gkliftzrho(jl,klev+1)2.romega
	1792	c * (2pstd(jl)+pmea(jl))
	1793	* 2pstd(jl)
	1794	* sin(rpi/180.plat(jl))pvlow(jl)
	1795	ztav(jl,klev+1)= gkliftzrho(jl,klev+1)2.romega
	1796	c * (2pstd(jl)+pmea(jl))
	1797	* 2pstd(jl)
	1798	* sin(rpi/180.plat(jl))pulow(jl)
	1799	else
	1800	ztau(jl,klev+1)=0.0
	1801	ztav(jl,klev+1)=0.0
	1802	endif
	1803	301 continue
	1804
	1805	C
	1806	C* 4. COMPUTE LIFT PROFILE
	1807	C* --------------------
	1808	C
	1809
	1810	400 continue
	1811
	1812	do 401 jk=1,klev
	1813	do 401 jl=kidia,kfdia
	1814	if(ktest(jl).eq.1) then
	1815	ztau(jl,jk)=ztau(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
	1816	ztav(jl,jk)=ztav(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
	1817	else
	1818	ztau(jl,jk)=0.0
	1819	ztav(jl,jk)=0.0
	1820	endif
	1821	401 continue
	1822	C
	1823	C
	1824	C* 5. COMPUTE TENDENCIES.
	1825	C* -------------------
	1826	if(lifthigh)then
	1827	C
	1828	500 continue
	1829	C
	1830	C EXPLICIT SOLUTION AT ALL LEVELS
	1831	C
	1832	do 524 jk=1,klev
	1833	do 523 jl=kidia,kfdia
	1834	if(ktest(jl).eq.1) then
	1835	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
	1836	zdudt(jl)=-rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
	1837	zdvdt(jl)=-rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
	1838	endif
	1839	523 continue
	1840	524 continue
	1841	C
	1842	C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
	1843	C
	1844	do 530 jk=1,klev
	1845	do 530 jl=kidia,kfdia
	1846	if(ktest(jl).eq.1) then
	1847
	1848	zslow=sqrt(pulow(jl)2+pvlow(jl)2)
	1849	zsqua=amax1(sqrt(pum1(jl,jk)2+pvm1(jl,jk)2),gvsec)
	1850	zscav=-zdudt(jl)pvm1(jl,jk)+zdvdt(jl)pum1(jl,jk)
	1851	if(zsqua.gt.gvsec)then
	1852	pvom(jl,jk)=-zscavpvm1(jl,jk)/zsqua*2
	1853	pvol(jl,jk)= zscavpum1(jl,jk)/zsqua*2
	1854	else
	1855	pvom(jl,jk)=0.0
	1856	pvol(jl,jk)=0.0
	1857	endif
	1858	zsqua=sqrt(pum1(jl,jk)2+pum1(jl,jk)2)
	1859	if(zsqua.lt.zslow)then
	1860	pvom(jl,jk)=zsqua/zslow*pvom(jl,jk)
	1861	pvol(jl,jk)=zsqua/zslow*pvol(jl,jk)
	1862	endif
	1863
	1864	endif
	1865	530 continue
	1866	C
	1867	C 6. LOW LEVEL LIFT, SEMI IMPLICIT:
	1868	C ----------------------------------
	1869
	1870	else
	1871
	1872	do 601 jl=kidia,kfdia
	1873	if(ktest(jl).eq.1) then
	1874	do jk=klev,iknub(jl),-1
	1875	zbet=gklift2.romegasin(rpi/180.plat(jl))ztmst
	1876	* (pgeom1(jl,iknub(jl)-1)-pgeom1(jl, jk))/
	1877	* (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
	1878	zdudt(jl)=-pum1(jl,jk)/ztmst/(1+zbet**2)
	1879	zdvdt(jl)=-pvm1(jl,jk)/ztmst/(1+zbet**2)
	1880	pvom(jl,jk)= zbet*2zdudt(jl) - zbet *zdvdt(jl)
	1881	pvol(jl,jk)= zbet zdudt(jl) + zbet2zdvdt(jl)
	1882	enddo
	1883	endif
	1884	601 continue
	1885
	1886	endif
	1887
	1888	c print *,' out orolift'
	1889
	1890	return
	1891	end
	1892	SUBROUTINE SUGWD_strato(NLON,NLEV,paprs,pplay)
	1893	C
	1894	C
	1895	C**** SUGWD INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
	1896	C
	1897	C PURPOSE.
	1898	C --------
	1899	C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
	1900	C GRAVITY WAVE DRAG PARAMETRIZATION.
	1901	C VERY IMPORTANT:
	1902	C ______________
	1903	C THIS ROUTINE SET_UP THE "TUNABLE PARAMETERS" OF THE
	1904	C VARIOUS SSO SCHEMES
	1905	C
	1906	C** INTERFACE.
	1907	C ----------
	1908	C CALL SUGWD FROM SUPHEC
	1909	C ----- ------
	1910	C
	1911	C EXPLICIT ARGUMENTS :
	1912	C --------------------
	1913	C PAPRS,PPLAY : Pressure at semi and full model levels
	1914	C NLEV : number of model levels
	1915	c NLON : number of points treated in the physics
	1916	C
	1917	C IMPLICIT ARGUMENTS :
	1918	C --------------------
	1919	C COMMON YOEGWD
	1920	C-GFRCRIT-R: Critical Non-dimensional mountain Height
	1921	C (HNC in (1), LOTT 1999)
	1922	C-GKWAKE--R: Bluff-body drag coefficient for low level wake
	1923	C (Cd in (2), LOTT 1999)
	1924	C-GRCRIT--R: Critical Richardson Number
	1925	C (Ric, End of first column p791 of LOTT 1999)
	1926	C-GKDRAG--R: Gravity wave drag coefficient
	1927	C (G in (3), LOTT 1999)
	1928	C-GKLIFT--R: Mountain Lift coefficient
	1929	C (Cl in (4), LOTT 1999)
	1930	C-GHMAX---R: Not used
	1931	C-GRAHILO-R: Set-up the trapped waves fraction
	1932	C (Beta , End of first column, LOTT 1999)
	1933	C
	1934	C-GSIGCR--R: Security value for blocked flow depth
	1935	C-NKTOPG--I: Security value for blocked flow level
	1936	C-nstra----I: An estimate to qualify the upper levels of
	1937	C the model where one wants to impose strees
	1938	C profiles
	1939	C-GSSECC--R: Security min value for low-level B-V frequency
	1940	C-GTSEC---R: Security min value for anisotropy and GW stress.
	1941	C-GVSEC---R: Security min value for ulow
	1942	C
	1943	C
	1944	C METHOD.
	1945	C -------
	1946	C SEE DOCUMENTATION
	1947	C
	1948	C EXTERNALS.
	1949	C ----------
	1950	C NONE
	1951	C
	1952	C REFERENCE.
	1953	C ----------
	1954	C Lott, 1999: Alleviation of stationary biases in a GCM through...
	1955	C Monthly Weather Review, 127, pp 788-801.
	1956	C
	1957	C AUTHOR.
	1958	C -------
	1959	C FRANCOIS LOTT LMD
	1960	C
	1961	C MODIFICATIONS.
	1962	C --------------
	1963	C ORIGINAL : 90-01-01 (MARTIN MILLER, ECMWF)
	1964	C LAST: 99-07-09 (FRANCOIS LOTT,LMD)
	1965	C ------------------------------------------------------------------
	1966	USE dimphy
	1967	USE mod_phys_lmdz_para
	1968	USE mod_grid_phy_lmdz
	1969	IMPLICIT NONE
	1970	C
	1971	C -----------------------------------------------------------------
	1972	#include "YOEGWD.h"
	1973	C ----------------------------------------------------------------
	1974	C
	1975	C ARGUMENTS
	1976	integer nlon,nlev
	1977	REAL paprs(nlon,nlev+1)
	1978	REAL pplay(nlon,nlev)
	1979	C
	1980	INTEGER JK
	1981	REAL ZPR,ZTOP,ZSIGT,ZPM1R
	1982	REAL :: pplay_glo(klon_glo,nlev)
	1983	REAL :: paprs_glo(klon_glo,nlev+1)
	1984
	1985	C
	1986	C* 1. SET THE VALUES OF THE PARAMETERS
	1987	C --------------------------------
	1988	C
	1989	100 CONTINUE
	1990	C
	1991	PRINT *,' DANS SUGWD NLEV=',NLEV
	1992	GHMAX=10000.
	1993	C
	1994	ZPR=100000.
	1995	ZTOP=0.001
	1996	ZSIGT=0.94
	1997	cold ZPR=80000.
	1998	cold ZSIGT=0.85
	1999	C
	2000	CALL gather(pplay,pplay_glo)
	2001	CALL bcast(pplay_glo)
	2002	CALL gather(paprs,paprs_glo)
	2003	CALL bcast(paprs_glo)
	2004
	2005	DO 110 JK=1,NLEV
[1491]	2006	ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2+1,1)
[1001]	2007	IF(ZPM1R.GE.ZSIGT)THEN
	2008	nktopg=JK
	2009	ENDIF
[1491]	2010	ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2+1,1)
[1001]	2011	IF(ZPM1R.GE.ZTOP)THEN
	2012	nstra=JK
	2013	ENDIF
	2014	110 CONTINUE
	2015	c
	2016	c inversion car dans orodrag on compte les niveaux a l'envers
	2017	nktopg=nlev-nktopg+1
	2018	nstra=nlev-nstra
	2019	print *,' DANS SUGWD nktopg=', nktopg
	2020	print *,' DANS SUGWD nstra=', nstra
	2021	C
	2022	GSIGCR=0.80
	2023	C
	2024	GKDRAG=0.1875
	2025	GRAHILO=0.1
	2026	GRCRIT=1.00
	2027	GFRCRIT=1.00
	2028	GKWAKE=0.50
	2029	C
	2030	GKLIFT=0.25
	2031	GVCRIT =0.1
	2032
	2033	WRITE(UNIT=6,FMT='('' * SSO essential constants *'')')
	2034	WRITE(UNIT=6,FMT='('' * SPECIFIED IN SUGWD *'')')
	2035	WRITE(UNIT=6,FMT='('' Gravity wave ct '',E13.7,'' '')')GKDRAG
	2036	WRITE(UNIT=6,FMT='('' Trapped/total wave dag '',E13.7,'' '')')
	2037	S GRAHILO
	2038	WRITE(UNIT=6,FMT='('' Critical Richardson = '',E13.7,'' '')')
	2039	S GRCRIT
	2040	WRITE(UNIT=6,FMT='('' Critical Froude'',e13.7)') GFRCRIT
	2041	WRITE(UNIT=6,FMT='('' Low level Wake bluff cte'',e13.7)') GKWAKE
	2042	WRITE(UNIT=6,FMT='('' Low level lift cte'',e13.7)') GKLIFT
	2043
	2044	C
	2045	C
	2046	C ----------------------------------------------------------------
	2047	C
	2048	C* 2. SET VALUES OF SECURITY PARAMETERS
	2049	C ---------------------------------
	2050	C
	2051	200 CONTINUE
	2052	C
	2053	GVSEC=0.10
	2054	GSSEC=0.0001
	2055	C
	2056	GTSEC=0.00001
	2057	C
	2058	RETURN
	2059	END
	2060

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

Download in other formats:

Original Format