Context Navigation

orografi_strato.F @ 3655

Last change on this file since 3655 was 1001, checked in by Laurent Fairhead, 16 years ago

Modifs sur le parallelisme: masquage dans la physique
Inclusion strato
mise en coherence etat0
le mode offline fonctionne maintenant en parallele,
les fichiers de la dynamiques sont correctement sortis et peuvent etre reconstruit avec rebuild
la version parallele de la dynamique peut s'executer sans MPI (sur 1 proc)
L'OPENMP fonctionne maintenant sans la parallelisation MPI.

YM
LF

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

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

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

Download in other formats:

Original Format