Context Navigation

orografi_strato.F @ 2206

Last change on this file since 2206 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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

File size: 57.6 KB

Line
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)
91	CHARACTER (LEN=20) :: modname='orografi_strato'
92	CHARACTER (LEN=80) :: abort_message
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
1684	CHARACTER (LEN=20) :: modname='orolift_strato'
1685	CHARACTER (LEN=80) :: abort_message
1686
1687
1688	C-----------------------------------------------------------------------
1689	C
1690	C* 1.1 initialisations
1691	C ---------------
1692
1693	lifthigh=.false.
1694
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
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
2006	ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2,1)
2007	IF(ZPM1R.GE.ZSIGT)THEN
2008	nktopg=JK
2009	ENDIF
2010	ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2,1)
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