Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

orografi_strato_mod.F90 @ 1

Last change on this file since 1 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.9 KB

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

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

Download in other formats:

Original Format