Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

orografi_strato.F90 @ 1939

Last change on this file since 1939 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 66.6 KB

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

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

Download in other formats:

Original Format