Context Navigation

orografi.F @ 5456

Last change on this file since 5456 was 634, checked in by Laurent Fairhead, 20 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 47.6 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE drag_noro (nlon,nlev,dtime,paprs,pplay,
5	e pmea,pstd, psig, pgam, pthe,ppic,pval,
6	e kgwd,kdx,ktest,
7	e t, u, v,
8	s pulow, pvlow, pustr, pvstr,
9	s d_t, d_u, d_v)
10	c
11	USE dimphy
12	IMPLICIT none
13	c======================================================================
14	c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
15	c Objet: Frottement de la montagne Interface
16	c======================================================================
17	c Arguments:
18	c dtime---input-R- pas d'integration (s)
19	c paprs---input-R-pression pour chaque inter-couche (en Pa)
20	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
21	c t-------input-R-temperature (K)
22	c u-------input-R-vitesse horizontale (m/s)
23	c v-------input-R-vitesse horizontale (m/s)
24	c
25	c d_t-----output-R-increment de la temperature
26	c d_u-----output-R-increment de la vitesse u
27	c d_v-----output-R-increment de la vitesse v
28	c======================================================================
29	cym#include "dimensions.h"
30	cym#include "dimphy.h"
31	#include "YOMCST.h"
32	c
33	c ARGUMENTS
34	c
35	INTEGER nlon,nlev
36	REAL dtime
37	REAL paprs(klon,klev+1)
38	REAL pplay(klon,klev)
39	REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
40	REAL ppic(nlon),pval(nlon)
41	REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
42	REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
43	REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
44	c
45	INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
46	c
47	c Variables locales:
48	c
49	REAL zgeom(klon,klev)
50	REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
51	REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
52	REAL papmf(klon,klev),papmh(klon,klev+1)
53	c
54	c initialiser les variables de sortie (pour securite)
55	c
56	DO i = 1,klon
57	pulow(i) = 0.0
58	pvlow(i) = 0.0
59	pustr(i) = 0.0
60	pvstr(i) = 0.0
61	ENDDO
62	DO k = 1, klev
63	DO i = 1, klon
64	d_t(i,k) = 0.0
65	d_u(i,k) = 0.0
66	d_v(i,k) = 0.0
67	pdudt(i,k)=0.0
68	pdvdt(i,k)=0.0
69	pdtdt(i,k)=0.0
70	ENDDO
71	ENDDO
72	c
73	c preparer les variables d'entree (attention: l'ordre des niveaux
74	c verticaux augmente du haut vers le bas)
75	c
76	DO k = 1, klev
77	DO i = 1, klon
78	pt(i,k) = t(i,klev-k+1)
79	pu(i,k) = u(i,klev-k+1)
80	pv(i,k) = v(i,klev-k+1)
81	papmf(i,k) = pplay(i,klev-k+1)
82	ENDDO
83	ENDDO
84	DO k = 1, klev+1
85	DO i = 1, klon
86	papmh(i,k) = paprs(i,klev-k+2)
87	ENDDO
88	ENDDO
89	DO i = 1, klon
90	zgeom(i,klev) = RD * pt(i,klev)
91	. * LOG(papmh(i,klev+1)/papmf(i,klev))
92	ENDDO
93	DO k = klev-1, 1, -1
94	DO i = 1, klon
95	zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
96	. * LOG(papmf(i,k+1)/papmf(i,k))
97	ENDDO
98	ENDDO
99	c
100	c appeler la routine principale
101	c
102	CALL orodrag(klon,klev,kgwd,kdx,ktest,
103	. dtime,
104	. papmh, papmf, zgeom,
105	. pt, pu, pv,
106	. pmea, pstd, psig, pgam, pthe, ppic,pval,
107	. pulow,pvlow,
108	. pdudt,pdvdt,pdtdt)
109	C
110	DO k = 1, klev
111	DO i = 1, klon
112	d_u(i,klev+1-k) = dtime*pdudt(i,k)
113	d_v(i,klev+1-k) = dtime*pdvdt(i,k)
114	d_t(i,klev+1-k) = dtime*pdtdt(i,k)
115	pustr(i) = pustr(i)
116	. +rgpdudt(i,k)(papmh(i,k+1)-papmh(i,k))
117	pvstr(i) = pvstr(i)
118	. +rgpdvdt(i,k)(papmh(i,k+1)-papmh(i,k))
119	ENDDO
120	ENDDO
121	c
122	RETURN
123	END
124	SUBROUTINE orodrag( nlon,nlev
125	i , kgwd, kdx, ktest
126	r , ptsphy
127	r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
128	r , pmea, pstd, psig, pgamma, ptheta, ppic, pval
129	c outputs
130	r , pulow,pvlow
131	r , pvom,pvol,pte )
132
133	USE dimphy
134	implicit none
135
136	c
137	c
138	c**** gwdrag - does the gravity wave parametrization.
139	c
140	c purpose.
141	c --------
142	c
143	c this routine computes the physical tendencies of the
144	c prognostic variables u,v and t due to vertical transports by
145	c subgridscale orographically excited gravity waves
146	c
147	c** interface.
148	c ----------
149	c called from callpar.
150	c
151	c the routine takes its input from the long-term storage:
152	c u,v,t and p at t-1.
153	c
154	c explicit arguments :
155	c --------------------
156	c ==== inputs ===
157	c ==== outputs ===
158	c
159	c implicit arguments : none
160	c --------------------
161	c
162	c implicit logical (l)
163	c
164	c method.
165	c -------
166	c
167	c externals.
168	c ----------
169	integer ismin, ismax
170	external ismin, ismax
171	c
172	c reference.
173	c ----------
174	c
175	c author.
176	c -------
177	c m.miller + b.ritter e.c.m.w.f. 15/06/86.
178	c
179	c f.lott + m. miller e.c.m.w.f. 22/11/94
180	c-----------------------------------------------------------------------
181	c
182	c
183	cym#include "dimensions.h"
184	cym#include "dimphy.h"
185	#include "YOMCST.h"
186	#include "YOEGWD.h"
187	c-----------------------------------------------------------------------
188	c
189	c* 0.1 arguments
190	c ---------
191	c
192	c
193	cym integer nlon, nlev, klevm1
194	integer nlon, nlev
195	integer kgwd, jl, ilevp1, jk, ji
196	real zdelp, ztemp, zforc, ztend
197	real rover, zb, zc, zconb, zabsv
198	real zzd1, ratio, zbet, zust,zvst, zdis
199	real pte(nlon,nlev),
200	* pvol(nlon,nlev),
201	* pvom(nlon,nlev),
202	* pulow(klon),
203	* pvlow(klon)
204	real pum1(nlon,nlev),
205	* pvm1(nlon,nlev),
206	* ptm1(nlon,nlev),
207	* pmea(nlon),pstd(nlon),psig(nlon),
208	* pgamma(nlon),ptheta(nlon),ppic(nlon),pval(nlon),
209	* pgeom1(nlon,nlev),
210	* papm1(nlon,nlev),
211	* paphm1(nlon,nlev+1)
212	c
213	integer kdx(nlon),ktest(nlon)
214	c-----------------------------------------------------------------------
215	c
216	c* 0.2 local arrays
217	c ------------
218	integer isect(klon),
219	* icrit(klon),
220	* ikcrith(klon),
221	* ikenvh(klon),
222	* iknu(klon),
223	* iknu2(klon),
224	* ikcrit(klon),
225	* ikhlim(klon)
226	c
227	real ztau(klon,klev+1),
228	$ ztauf(klon,klev+1),
229	* zstab(klon,klev+1),
230	* zvph(klon,klev+1),
231	* zrho(klon,klev+1),
232	* zri(klon,klev+1),
233	* zpsi(klon,klev+1),
234	* zzdep(klon,klev)
235	real zdudt(klon),
236	* zdvdt(klon),
237	* zdtdt(klon),
238	* zdedt(klon),
239	* zvidis(klon),
240	* znu(klon),
241	* zd1(klon),
242	* zd2(klon),
243	* zdmod(klon)
244	real ztmst, ptsphy, zrtmst
245	c
246	c------------------------------------------------------------------
247	c
248	c* 1. initialization
249	c --------------
250	c
251	100 continue
252	c
253	c ------------------------------------------------------------------
254	c
255	c* 1.1 computational constants
256	c -----------------------
257	c
258	110 continue
259	c
260	c ztmst=twodt
261	c if(nstep.eq.nstart) ztmst=0.5*twodt
262	cym klevm1=klev-1
263	ztmst=ptsphy
264	zrtmst=1./ztmst
265	c ------------------------------------------------------------------
266	c
267	120 continue
268	c
269	c ------------------------------------------------------------------
270	c
271	c* 1.3 check whether row contains point for printing
272	c ---------------------------------------------
273	c
274	130 continue
275	c
276	c ------------------------------------------------------------------
277	c
278	c* 2. precompute basic state variables.
279	c* ---------- ----- ----- ----------
280	c* define low level wind, project winds in plane of
281	c* low level wind, determine sector in which to take
282	c* the variance and set indicator for critical levels.
283	c
284	200 continue
285	c
286	c
287	c
288	call orosetup
289	* ( nlon, ktest
290	* , ikcrit, ikcrith, icrit, ikenvh,iknu,iknu2
291	* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
292	* , zrho , zri , zstab , ztau , zvph , zpsi, zzdep
293	* , pulow, pvlow
294	* , ptheta,pgamma,pmea,ppic,pval,znu ,zd1, zd2, zdmod )
295	c
296	c
297	c
298	c***********************************************************
299	c
300	c
301	c* 3. compute low level stresses using subcritical and
302	c* supercritical forms.computes anisotropy coefficient
303	c* as measure of orographic twodimensionality.
304	c
305	300 continue
306	c
307	call gwstress
308	* ( nlon , nlev
309	* , ktest , icrit, ikenvh, iknu
310	* , zrho , zstab, zvph , pstd, psig, pmea, ppic
311	* , ztau
312	* , pgeom1,zdmod)
313	c
314	c
315	c* 4. compute stress profile.
316	c* ------- ------ --------
317	c
318	400 continue
319	c
320	c
321	call gwprofil
322	* ( nlon , nlev
323	* , kgwd , kdx , ktest
324	* , ikcrith, icrit
325	* , paphm1, zrho , zstab , zvph
326	* , zri , ztau
327	* , zdmod , psig , pstd)
328	c
329	c
330	c* 5. compute tendencies.
331	c* -------------------
332	c
333	500 continue
334	c
335	c explicit solution at all levels for the gravity wave
336	c implicit solution for the blocked levels
337
338	do 510 jl=kidia,kfdia
339	zvidis(jl)=0.0
340	zdudt(jl)=0.0
341	zdvdt(jl)=0.0
342	zdtdt(jl)=0.0
343	510 continue
344	c
345	ilevp1=klev+1
346	c
347	c
348	do 524 jk=1,klev
349	c
350	c
351	c do 523 jl=1,kgwd
352	c ji=kdx(jl)
353	c Modif vectorisation 02/04/2004
354	do 523 ji=kidia,kfdia
355	if(ktest(ji).eq.1) then
356
357	zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
358	ztemp=-rg(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)zdelp)
359	zdudt(ji)=(pulow(ji)zd1(ji)-pvlow(ji)zd2(ji))*ztemp/zdmod(ji)
360	zdvdt(ji)=(pvlow(ji)zd1(ji)+pulow(ji)zd2(ji))*ztemp/zdmod(ji)
361	c
362	c controle des overshoots:
363	c
364	zforc=sqrt(zdudt(ji)2+zdvdt(ji)2)+1.E-12
365	ztend=sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/ztmst+1.E-12
366	rover=0.25
367	if(zforc.ge.rover*ztend)then
368	zdudt(ji)=roverztend/zforczdudt(ji)
369	zdvdt(ji)=roverztend/zforczdvdt(ji)
370	endif
371	c
372	c fin du controle des overshoots
373	c
374	if(jk.ge.ikenvh(ji)) then
375	zb=1.0-0.18pgamma(ji)-0.04pgamma(ji)**2
376	zc=0.48pgamma(ji)+0.3pgamma(ji)**2
377	zconb=2.ztmstgkwakepsig(ji)/(4.pstd(ji))
378	zabsv=sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/2.
379	zzd1=zbcos(zpsi(ji,jk))2+zcsin(zpsi(ji,jk))**2
380	ratio=(cos(zpsi(ji,jk))*2+pgamma(ji)sin(zpsi(ji,jk))**2)/
381	* (pgamma(ji)cos(zpsi(ji,jk))2+sin(zpsi(ji,jk))*2)
382	zbet=max(0.,2.-1./ratio)zconbzzdep(ji,jk)zzd1zabsv
383	c
384	c simplement oppose au vent
385	c
386	zdudt(ji)=-pum1(ji,jk)/ztmst
387	zdvdt(ji)=-pvm1(ji,jk)/ztmst
388	c
389	c projection dans la direction de l'axe principal de l'orographie
390	cmod zdudt(ji)=-(pum1(ji,jk)cos(ptheta(ji)rpi/180.)
391	cmod * +pvm1(ji,jk)sin(ptheta(ji)rpi/180.))
392	cmod * cos(ptheta(ji)rpi/180.)/ztmst
393	cmod zdvdt(ji)=-(pum1(ji,jk)cos(ptheta(ji)rpi/180.)
394	cmod * +pvm1(ji,jk)sin(ptheta(ji)rpi/180.))
395	cmod * sin(ptheta(ji)rpi/180.)/ztmst
396	zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
397	zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
398	end if
399	pvom(ji,jk)=zdudt(ji)
400	pvol(ji,jk)=zdvdt(ji)
401	zust=pum1(ji,jk)+ztmst*zdudt(ji)
402	zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
403	zdis=0.5(pum1(ji,jk)2+pvm1(ji,jk)2-zust2-zvst*2)
404	zdedt(ji)=zdis/ztmst
405	zvidis(ji)=zvidis(ji)+zdis*zdelp
406	zdtdt(ji)=zdedt(ji)/rcpd
407	c pte(ji,jk)=zdtdt(ji)
408	c
409	c ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
410	c
411	pte(ji,jk)=0.0
412
413	endif
414	523 continue
415
416	524 continue
417	c
418	c
419	return
420	end
421	SUBROUTINE orosetup
422	* ( nlon , ktest
423	* , kkcrit, kkcrith, kcrit
424	* , kkenvh, kknu , kknu2
425	* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
426	* , prho , pri , pstab , ptau , pvph ,ppsi, pzdep
427	* , pulow , pvlow
428	* , ptheta, pgamma, pmea, ppic, pval
429	* , pnu , pd1 , pd2 ,pdmod )
430	c
431	c**** gwsetup
432	c
433	c purpose.
434	c --------
435	c
436	c** interface.
437	c ----------
438	c from orodrag
439	c
440	c explicit arguments :
441	c --------------------
442	c ==== inputs ===
443	c ==== outputs ===
444	c
445	c implicit arguments : none
446	c --------------------
447	c
448	c method.
449	c -------
450	c
451	c
452	c externals.
453	c ----------
454	c
455	c
456	c reference.
457	c ----------
458	c
459	c see ecmwf research department documentation of the "i.f.s."
460	c
461	c author.
462	c -------
463	c
464	c modifications.
465	c --------------
466	c f.lott for the new-gwdrag scheme november 1993
467	c
468	c-----------------------------------------------------------------------
469	USE dimphy
470	implicit none
471	c
472
473	cym#include "dimensions.h"
474	cym#include "dimphy.h"
475	#include "YOMCST.h"
476	#include "YOEGWD.h"
477
478	c-----------------------------------------------------------------------
479	c
480	c* 0.1 arguments
481	c ---------
482	c
483	integer nlon
484	integer jl, jk
485	real zdelp
486
487	integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),
488	* ktest(nlon),kkenvh(nlon)
489
490	c
491	real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
492	* pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
493	* prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
494	* ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
495	* pzdep(nlon,klev)
496	real pulow(nlon),pvlow(nlon),ptheta(nlon),pgamma(nlon),pnu(nlon),
497	* pd1(nlon),pd2(nlon),pdmod(nlon)
498	real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
499	c
500	c-----------------------------------------------------------------------
501	c
502	c* 0.2 local arrays
503	c ------------
504	c
505	c
506	integer ilevm1, ilevm2, ilevh
507	real zcons1, zcons2,zcons3, zhgeo
508	real zu, zphi, zvt1,zvt2, zst, zvar, zdwind, zwind
509	real zstabm, zstabp, zrhom, zrhop, alpha
510	real zggeenv, zggeom1,zgvar
511	logical lo
512	logical ll1(klon,klev+1)
513	integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
514	* kentp(klon),ncount(klon)
515	c
516	real zhcrit(klon,klev),zvpf(klon,klev),
517	* zdp(klon,klev)
518	real znorm(klon),zb(klon),zc(klon),
519	* zulow(klon),zvlow(klon),znup(klon),znum(klon)
520	c
521	c ------------------------------------------------------------------
522	c
523	c* 1. initialization
524	c --------------
525	c
526	c print *,' entree gwsetup'
527	100 continue
528	c
529	c ------------------------------------------------------------------
530	c
531	c* 1.1 computational constants
532	c -----------------------
533	c
534	110 continue
535	c
536	ilevm1=klev-1
537	ilevm2=klev-2
538	ilevh =klev/3
539	c
540	zcons1=1./rd
541	cold zcons2=g**2/cpd
542	zcons2=rg**2/rcpd
543	cold zcons3=1.5*api
544	zcons3=1.5*rpi
545	c
546	c
547	c ------------------------------------------------------------------
548	c
549	c* 2.
550	c --------------
551	c
552	200 continue
553	c
554	c ------------------------------------------------------------------
555	c
556	c* 2.1 define low level wind, project winds in plane of
557	c* low level wind, determine sector in which to take
558	c* the variance and set indicator for critical levels.
559	c
560	c
561	c
562	do 2001 jl=kidia,kfdia
563	kknu(jl) =klev
564	kknu2(jl) =klev
565	kknub(jl) =klev
566	kknul(jl) =klev
567	pgamma(jl) =max(pgamma(jl),gtsec)
568	ll1(jl,klev+1)=.false.
569	2001 continue
570	c
571	c Ajouter une initialisation (L. Li, le 23fev99):
572	c
573	do jk=klev,ilevh,-1
574	do jl=kidia,kfdia
575	ll1(jl,jk)= .FALSE.
576	ENDDO
577	ENDDO
578	c
579	c* define top of low level flow
580	c ----------------------------
581	do 2002 jk=klev,ilevh,-1
582	do 2003 jl=kidia,kfdia
583	lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
584	if(lo) then
585	kkcrit(jl)=jk
586	endif
587	zhcrit(jl,jk)=ppic(jl)
588	zhgeo=pgeom1(jl,jk)/rg
589	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
590	if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
591	kknu(jl)=jk
592	endif
593	if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
594	2003 continue
595	2002 continue
596	do 2004 jk=klev,ilevh,-1
597	do 2005 jl=kidia,kfdia
598	zhcrit(jl,jk)=ppic(jl)-pval(jl)
599	zhgeo=pgeom1(jl,jk)/rg
600	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
601	if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
602	kknu2(jl)=jk
603	endif
604	if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
605	2005 continue
606	2004 continue
607	do 2006 jk=klev,ilevh,-1
608	do 2007 jl=kidia,kfdia
609	zhcrit(jl,jk)=amax1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
610	zhgeo=pgeom1(jl,jk)/rg
611	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
612	if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
613	kknub(jl)=jk
614	endif
615	if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
616	2007 continue
617	2006 continue
618	c
619	do 2010 jl=kidia,kfdia
620	kknu(jl)=min(kknu(jl),nktopg)
621	kknu2(jl)=min(kknu2(jl),nktopg)
622	kknub(jl)=min(kknub(jl),nktopg)
623	kknul(jl)=klev
624	2010 continue
625	c
626
627	210 continue
628	c
629	c
630	cc* initialize various arrays
631	c
632	do 2107 jl=kidia,kfdia
633	prho(jl,klev+1) =0.0
634	pstab(jl,klev+1) =0.0
635	pstab(jl,1) =0.0
636	pri(jl,klev+1) =9999.0
637	ppsi(jl,klev+1) =0.0
638	pri(jl,1) =0.0
639	pvph(jl,1) =0.0
640	pulow(jl) =0.0
641	pvlow(jl) =0.0
642	zulow(jl) =0.0
643	zvlow(jl) =0.0
644	kkcrith(jl) =klev
645	kkenvh(jl) =klev
646	kentp(jl) =klev
647	kcrit(jl) =1
648	ncount(jl) =0
649	ll1(jl,klev+1) =.false.
650	2107 continue
651	c
652	c* define low-level flow
653	c ---------------------
654	c
655	do 223 jk=klev,2,-1
656	do 222 jl=kidia,kfdia
657	if(ktest(jl).eq.1) then
658	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
659	prho(jl,jk)=2.paphm1(jl,jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
660	pstab(jl,jk)=2.zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))
661	* (1.-rcpdprho(jl,jk)(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
662	pstab(jl,jk)=max(pstab(jl,jk),gssec)
663	endif
664	222 continue
665	223 continue
666	c
667	c********************************************************************
668	c
669	c* define blocked flow
670	c -------------------
671	do 2115 jk=klev,ilevh,-1
672	do 2116 jl=kidia,kfdia
673	if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
674	pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
675	pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
676	end if
677	2116 continue
678	2115 continue
679	do 2110 jl=kidia,kfdia
680	pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
681	pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
682	znorm(jl)=max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
683	pvph(jl,klev+1)=znorm(jl)
684	2110 continue
685	c
686	c***** setup orography axes and define plane of profiles *****
687	c
688	do 2112 jl=kidia,kfdia
689	lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
690	if(lo) then
691	zu=pulow(jl)+2.*gvsec
692	else
693	zu=pulow(jl)
694	endif
695	zphi=atan(pvlow(jl)/zu)
696	ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
697	zb(jl)=1.-0.18pgamma(jl)-0.04pgamma(jl)**2
698	zc(jl)=0.48pgamma(jl)+0.3pgamma(jl)**2
699	pd1(jl)=zb(jl)-(zb(jl)-zc(jl))(sin(ppsi(jl,klev+1))*2)
700	pd2(jl)=(zb(jl)-zc(jl))sin(ppsi(jl,klev+1))cos(ppsi(jl,klev+1))
701	pdmod(jl)=sqrt(pd1(jl)2+pd2(jl)2)
702	2112 continue
703	c
704	c ********** define flow in plane of lowlevel stress ***********
705	c
706	do 213 jk=1,klev
707	do 212 jl=kidia,kfdia
708	if(ktest(jl).eq.1) then
709	zvt1 =pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk)
710	zvt2 =-pvlow(jl)pum1(jl,jk)+pulow(jl)pvm1(jl,jk)
711	zvpf(jl,jk)=(zvt1pd1(jl)+zvt2pd2(jl))/(znorm(jl)*pdmod(jl))
712	endif
713	ptau(jl,jk) =0.0
714	pzdep(jl,jk) =0.0
715	ppsi(jl,jk) =0.0
716	ll1(jl,jk) =.false.
717	212 continue
718	213 continue
719	do 215 jk=2,klev
720	do 214 jl=kidia,kfdia
721	if(ktest(jl).eq.1) then
722	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
723	pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
724	* (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
725	* /zdp(jl,jk)
726	if(pvph(jl,jk).lt.gvsec) then
727	pvph(jl,jk)=gvsec
728	kcrit(jl)=jk
729	endif
730	endif
731	214 continue
732	215 continue
733	c
734	c
735	c* 2.2 brunt-vaisala frequency and density at half levels.
736	c
737	220 continue
738	c
739	do 2211 jk=ilevh,klev
740	do 221 jl=kidia,kfdia
741	if(ktest(jl).eq.1) then
742	if(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) then
743	zst=zcons2/ptm1(jl,jk)(1.-rcpdprho(jl,jk)*
744	* (ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
745	pstab(jl,klev+1)=pstab(jl,klev+1)+zst*zdp(jl,jk)
746	pstab(jl,klev+1)=max(pstab(jl,klev+1),gssec)
747	prho(jl,klev+1)=prho(jl,klev+1)+paphm1(jl,jk)2.zdp(jl,jk)
748	* *zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
749	endif
750	endif
751	221 continue
752	2211 continue
753	c
754	do 2212 jl=kidia,kfdia
755	pstab(jl,klev+1)=pstab(jl,klev+1)/(papm1(jl,kknul(jl))
756	* -papm1(jl,kknub(jl)))
757	prho(jl,klev+1)=prho(jl,klev+1)/(papm1(jl,kknul(jl))
758	* -papm1(jl,kknub(jl)))
759	zvar=pstd(jl)
760	2212 continue
761	c
762	c* 2.3 mean flow richardson number.
763	c* and critical height for froude layer
764	c
765	230 continue
766	c
767	do 232 jk=2,klev
768	do 231 jl=kidia,kfdia
769	if(ktest(jl).eq.1) then
770	zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
771	pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
772	* /(rgprho(jl,jk)zdwind))**2
773	pri(jl,jk)=max(pri(jl,jk),grcrit)
774	endif
775	231 continue
776	232 continue
777
778	c
779	c
780	c* define top of 'envelope' layer
781	c ----------------------------
782
783	do 233 jl=kidia,kfdia
784	pnu (jl)=0.0
785	znum(jl)=0.0
786	233 continue
787
788	do 234 jk=2,klev-1
789	do 234 jl=kidia,kfdia
790
791	if(ktest(jl).eq.1) then
792
793	if (jk.ge.kknub(jl)) then
794
795	znum(jl)=pnu(jl)
796	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
797	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
798	zwind=max(sqrt(zwind**2),gvsec)
799	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
800	zstabm=sqrt(max(pstab(jl,jk ),gssec))
801	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
802	zrhom=prho(jl,jk )
803	zrhop=prho(jl,jk+1)
804	pnu(jl) = pnu(jl) + (zdelp/rg)*
805	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
806	if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
807	* .and.(kkenvh(jl).eq.klev))
808	* kkenvh(jl)=jk
809
810	endif
811
812	endif
813
814	234 continue
815
816	c calculation of a dynamical mixing height for the breaking
817	c of gravity waves:
818
819
820	do 235 jl=kidia,kfdia
821	znup(jl)=0.0
822	znum(jl)=0.0
823	235 continue
824
825	do 236 jk=klev-1,2,-1
826	do 236 jl=kidia,kfdia
827
828	if(ktest(jl).eq.1) then
829
830	znum(jl)=znup(jl)
831	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
832	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
833	zwind=max(sqrt(zwind**2),gvsec)
834	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
835	zstabm=sqrt(max(pstab(jl,jk ),gssec))
836	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
837	zrhom=prho(jl,jk )
838	zrhop=prho(jl,jk+1)
839	znup(jl) = znup(jl) + (zdelp/rg)*
840	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
841	if((znum(jl).le.rpi/2.).and.(znup(jl).gt.rpi/2.)
842	* .and.(kkcrith(jl).eq.klev))
843	* kkcrith(jl)=jk
844
845	endif
846
847	236 continue
848
849	do 237 jl=kidia,kfdia
850	kkcrith(jl)=min0(kkcrith(jl),kknu2(jl))
851	kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
852	237 continue
853	c
854	c directional info for flow blocking *************************
855	c
856	do 251 jk=ilevh,klev
857	do 252 jl=kidia,kfdia
858	if(jk.ge.kkenvh(jl)) then
859	lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
860	if(lo) then
861	zu=pum1(jl,jk)+2.*gvsec
862	else
863	zu=pum1(jl,jk)
864	endif
865	zphi=atan(pvm1(jl,jk)/zu)
866	ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
867	end if
868	252 continue
869	251 continue
870	c forms the vertical 'leakiness' **************************
871
872	alpha=3.
873
874	do 254 jk=ilevh,klev
875	do 253 jl=kidia,kfdia
876	if(jk.ge.kkenvh(jl)) then
877	zggeenv=amax1(1.,
878	* (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)
879	zggeom1=amax1(pgeom1(jl,jk),1.)
880	zgvar=amax1(pstd(jl)*rg,1.)
881	cmod pzdep(jl,jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar))
882	pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl, jk))/
883	* (pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,klev))
884	end if
885	253 continue
886	254 continue
887
888	260 continue
889
890	return
891	end
892	SUBROUTINE gwstress
893	* ( nlon , nlev
894	* , ktest, kcrit, kkenvh
895	* , kknu
896	* , prho , pstab , pvph , pstd, psig
897	* , pmea , ppic , ptau
898	* , pgeom1 , pdmod )
899	c
900	c**** gwstress
901	c
902	c purpose.
903	c --------
904	c
905	c** interface.
906	c ----------
907	c call gwstress from gwdrag
908	c
909	c explicit arguments :
910	c --------------------
911	c ==== inputs ===
912	c ==== outputs ===
913	c
914	c implicit arguments : none
915	c --------------------
916	c
917	c method.
918	c -------
919	c
920	c
921	c externals.
922	c ----------
923	c
924	c
925	c reference.
926	c ----------
927	c
928	c see ecmwf research department documentation of the "i.f.s."
929	c
930	c author.
931	c -------
932	c
933	c modifications.
934	c --------------
935	c f. lott put the new gwd on ifs 22/11/93
936	c
937	c-----------------------------------------------------------------------
938	USE dimphy
939	implicit none
940	cym#include "dimensions.h"
941	cym#include "dimphy.h"
942	#include "YOMCST.h"
943	#include "YOEGWD.h"
944
945	c-----------------------------------------------------------------------
946	c
947	c* 0.1 arguments
948	c ---------
949	c
950	integer nlon, nlev
951	integer kcrit(nlon),
952	* ktest(nlon),kkenvh(nlon),kknu(nlon)
953	c
954	real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1),
955	* pvph(nlon,nlev+1),
956	* pgeom1(nlon,nlev),pstd(nlon)
957	c
958	real psig(nlon)
959	real pmea(nlon),ppic(nlon)
960	real pdmod(nlon)
961	c
962	c-----------------------------------------------------------------------
963	c
964	c* 0.2 local arrays
965	c ------------
966	integer jl
967	real zblock, zvar, zeff
968	logical lo
969	c
970	c-----------------------------------------------------------------------
971	c
972	c* 0.3 functions
973	c ---------
974	c ------------------------------------------------------------------
975	c
976	c* 1. initialization
977	c --------------
978	c
979	100 continue
980	c
981	c* 3.1 gravity wave stress.
982	c
983	300 continue
984	c
985	c
986	do 301 jl=kidia,kfdia
987	if(ktest(jl).eq.1) then
988
989	c effective mountain height above the blocked flow
990
991	if(kkenvh(jl).eq.klev)then
992	zblock=0.0
993	else
994	zblock=(pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg
995	endif
996
997	zvar=ppic(jl)-pmea(jl)
998	zeff=amax1(0.,zvar-zblock)
999
1000	ptau(jl,klev+1)=prho(jl,klev+1)gkdragpsig(jl)zeff*2
1001	* /4./pstd(jl)pvph(jl,klev+1)pdmod(jl)*sqrt(pstab(jl,klev+1))
1002
1003	c too small value of stress or low level flow include critical level
1004	c or low level flow: gravity wave stress nul.
1005
1006	lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl))
1007	* .or.(pvph(jl,klev+1).lt.gvcrit)
1008	c if(lo) ptau(jl,klev+1)=0.0
1009
1010	else
1011
1012	ptau(jl,klev+1)=0.0
1013
1014	endif
1015
1016	301 continue
1017	c
1018	return
1019	end
1020	SUBROUTINE GWPROFIL
1021	* ( NLON, NLEV
1022	* , kgwd, kdx , ktest
1023	* , KKCRITH, KCRIT
1024	* , PAPHM1, PRHO , PSTAB , PVPH , PRI , PTAU
1025	* , pdmod , psig , pvar)
1026
1027	C**** GWPROFIL
1028	C
1029	C PURPOSE.
1030	C --------
1031	C
1032	C** INTERFACE.
1033	C ----------
1034	C FROM GWDRAG
1035	C
1036	C EXPLICIT ARGUMENTS :
1037	C --------------------
1038	C ==== INPUTS ===
1039	C ==== OUTPUTS ===
1040	C
1041	C IMPLICIT ARGUMENTS : NONE
1042	C --------------------
1043	C
1044	C METHOD:
1045	C -------
1046	C THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
1047	C IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
1048	C AND THE TOP OF THE BLOCKED LAYER (KKENVH).
1049	C IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE
1050	C BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR
1051	C NONLINEAR GRAVITY WAVE BREAKING.
1052	C ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
1053	C LEVEL (KCRIT) OR WHEN IT BREAKS.
1054	C
1055	C
1056	C
1057	C EXTERNALS.
1058	C ----------
1059	C
1060	C
1061	C REFERENCE.
1062	C ----------
1063	C
1064	C SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
1065	C
1066	C AUTHOR.
1067	C -------
1068	C
1069	C MODIFICATIONS.
1070	C --------------
1071	C PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
1072	C-----------------------------------------------------------------------
1073	USE dimphy
1074	implicit none
1075	C
1076
1077	C
1078
1079	cym#include "dimensions.h"
1080	cym#include "dimphy.h"
1081	#include "YOMCST.h"
1082	#include "YOEGWD.h"
1083
1084	C-----------------------------------------------------------------------
1085	C
1086	C* 0.1 ARGUMENTS
1087	C ---------
1088	C
1089	integer nlon,nlev
1090	INTEGER KKCRITH(NLON),KCRIT(NLON)
1091	* ,kdx(nlon) , ktest(nlon)
1092
1093	C
1094	REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1),
1095	* PRHO (NLON,NLEV+1), PVPH (NLON,NLEV+1),
1096	* PRI (NLON,NLEV+1), PTAU(NLON,NLEV+1)
1097
1098	REAL pdmod (NLON) , psig(NLON),
1099	* pvar(NLON)
1100
1101	C-----------------------------------------------------------------------
1102	C
1103	C* 0.2 LOCAL ARRAYS
1104	C ------------
1105	C
1106	integer ilevh, ji, kgwd, jl, jk
1107	real zsqr, zalfa, zriw, zdel, zb, zalpha,zdz2n
1108	real zdelp, zdelpt
1109	REAL ZDZ2 (KLON,KLEV) , ZNORM(KLON) , zoro(KLON)
1110	REAL ZTAU (KLON,KLEV+1)
1111	C
1112	C-----------------------------------------------------------------------
1113	C
1114	C* 1. INITIALIZATION
1115	C --------------
1116	C
1117	c print *,' entree gwprofil'
1118	100 CONTINUE
1119	C
1120	C
1121	C* COMPUTATIONAL CONSTANTS.
1122	C ------------- ----------
1123	C
1124	ilevh=KLEV/3
1125	C
1126	c DO 400 ji=1,kgwd
1127	c jl=kdx(ji)
1128	c Modif vectorisation 02/04/2004
1129	DO 400 jl=kidia,kfdia
1130	if (ktest(jl).eq.1) then
1131	Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0)
1132	ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1)
1133	endif
1134	400 CONTINUE
1135
1136	C
1137	DO 430 JK=KLEV,2,-1
1138	C
1139	C
1140	C* 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
1141	C BLOCKING LAYER.
1142	410 CONTINUE
1143	C
1144	c DO 411 ji=1,kgwd
1145	c jl=kdx(ji)
1146	c Modif vectorisation 02/04/2004
1147	do 411 jl=kidia,kfdia
1148	if (ktest(jl).eq.1) then
1149	IF(JK.GT.KKCRITH(JL)) THEN
1150	PTAU(JL,JK)=ZTAU(JL,KLEV+1)
1151	C ENDIF
1152	C IF(JK.EQ.KKCRITH(JL)) THEN
1153	ELSE
1154	PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1)
1155	ENDIF
1156	endif
1157	411 CONTINUE
1158	C
1159	C* 4.15 CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
1160	C LOW LEVEL FLOW LAYER.
1161	415 CONTINUE
1162	C
1163	C
1164	C* 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
1165	C
1166	420 CONTINUE
1167	C
1168	c DO 421 ji=1,kgwd
1169	c jl=kdx(ji)
1170	c Modif vectorisation 02/04/2004
1171	do 421 jl=kidia,kfdia
1172	if(ktest(jl).eq.1) then
1173	IF(JK.LT.KKCRITH(JL)) THEN
1174	ZNORM(JL)=gkdragPRHO(JL,JK)SQRT(PSTAB(JL,JK))*PVPH(JL,JK)
1175	* *zoro(jl)
1176	ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec)
1177	ENDIF
1178	endif
1179	421 CONTINUE
1180	C
1181	C* 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
1182	C* AND STRESS: BREAKING EVALUATION AND CRITICAL
1183	C LEVEL
1184	C
1185
1186	c DO 431 ji=1,kgwd
1187	c jl=Kdx(ji)
1188	c Modif vectorisation 02/04/2004
1189	do 431 jl=kidia,kfdia
1190	if(ktest(jl).eq.1) then
1191
1192	IF(JK.LT.KKCRITH(JL)) THEN
1193	IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN
1194	PTAU(JL,JK)=0.0
1195	ELSE
1196	ZSQR=SQRT(PRI(JL,JK))
1197	ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK)
1198	ZRIW=PRI(JL,JK)(1.-ZALFA)/(1+ZALFAZSQR)**2
1199	IF(ZRIW.LT.GRCRIT) THEN
1200	ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
1201	ZB=1./GRCRIT+2./ZSQR
1202	ZALPHA=0.5*(-ZB+SQRT(ZDEL))
1203	ZDZ2N=(PVPH(JL,JK)ZALPHA)*2/PSTAB(JL,JK)
1204	PTAU(JL,JK)=ZNORM(JL)*ZDZ2N
1205	ELSE
1206	PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
1207	ENDIF
1208	PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1))
1209	ENDIF
1210	ENDIF
1211	endif
1212	431 CONTINUE
1213
1214	430 CONTINUE
1215	440 CONTINUE
1216
1217	C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
1218
1219	c DO 530 ji=1,kgwd
1220	c jl=kdx(ji)
1221	c Modif vectorisation 02/04/2004
1222	do 530 jl=kidia,kfdia
1223	if(ktest(jl).eq.1) then
1224	ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL))
1225	ZTAU(JL,NSTRA)=PTAU(JL,NSTRA)
1226	endif
1227	530 CONTINUE
1228
1229	DO 531 JK=1,KLEV
1230
1231	c DO 532 ji=1,kgwd
1232	c jl=kdx(ji)
1233	c Modif vectorisation 02/04/2004
1234	do 532 jl=kidia,kfdia
1235	if(ktest(jl).eq.1) then
1236
1237
1238	IF(JK.GT.KKCRITH(JL))THEN
1239
1240	ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KLEV+1 )
1241	ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KLEV+1 )
1242	PTAU(JL,JK)=ZTAU(JL,KLEV+1 ) +
1243	. (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KLEV+1 ) )*
1244	. ZDELP/ZDELPT
1245
1246	ENDIF
1247
1248	endif
1249	532 CONTINUE
1250
1251	C REORGANISATION IN THE STRATOSPHERE
1252
1253	c DO 533 ji=1,kgwd
1254	c jl=kdx(ji)
1255	c Modif vectorisation 02/04/2004
1256	do 533 jl=kidia,kfdia
1257	if(ktest(jl).eq.1) then
1258
1259
1260	IF(JK.LT.NSTRA)THEN
1261
1262	ZDELP =PAPHM1(JL,NSTRA)
1263	ZDELPT=PAPHM1(JL,JK)
1264	PTAU(JL,JK)=ZTAU(JL,NSTRA)*ZDELPT/ZDELP
1265
1266	ENDIF
1267
1268	endif
1269	533 CONTINUE
1270
1271	C REORGANISATION IN THE TROPOSPHERE
1272
1273	c DO 534 ji=1,kgwd
1274	c jl=kdx(ji)
1275	c Modif vectorisation 02/04/2004
1276	do 534 jl=kidia,kfdia
1277	if(ktest(jl).eq.1) then
1278
1279
1280	IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN
1281
1282	ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKCRITH(JL))
1283	ZDELPT=PAPHM1(JL,NSTRA)-PAPHM1(JL,KKCRITH(JL))
1284	PTAU(JL,JK)=ZTAU(JL,KKCRITH(JL)) +
1285	* (ZTAU(JL,NSTRA)-ZTAU(JL,KKCRITH(JL)))*ZDELP
1286	. /ZDELPT
1287
1288	ENDIF
1289	endif
1290	534 CONTINUE
1291
1292
1293	531 CONTINUE
1294
1295
1296	RETURN
1297	END
1298	SUBROUTINE lift_noro (nlon,nlev,dtime,paprs,pplay,
1299	e plat,pmea,pstd, ppic,
1300	e ktest,
1301	e t, u, v,
1302	s pulow, pvlow, pustr, pvstr,
1303	s d_t, d_u, d_v)
1304	c
1305	USE dimphy
1306	IMPLICIT none
1307	c======================================================================
1308	c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
1309	c Objet: Frottement de la montagne Interface
1310	c======================================================================
1311	c Arguments:
1312	c dtime---input-R- pas d'integration (s)
1313	c paprs---input-R-pression pour chaque inter-couche (en Pa)
1314	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
1315	c t-------input-R-temperature (K)
1316	c u-------input-R-vitesse horizontale (m/s)
1317	c v-------input-R-vitesse horizontale (m/s)
1318	c
1319	c d_t-----output-R-increment de la temperature
1320	c d_u-----output-R-increment de la vitesse u
1321	c d_v-----output-R-increment de la vitesse v
1322	c======================================================================
1323	cym#include "dimensions.h"
1324	cym#include "dimphy.h"
1325	#include "YOMCST.h"
1326	c
1327	c ARGUMENTS
1328	c
1329	INTEGER nlon,nlev
1330	REAL dtime
1331	REAL paprs(klon,klev+1)
1332	REAL pplay(klon,klev)
1333	REAL plat(nlon),pmea(nlon)
1334	REAL pstd(nlon)
1335	REAL ppic(nlon)
1336	REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
1337	REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
1338	REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
1339	c
1340	INTEGER i, k, ktest(nlon)
1341	c
1342	c Variables locales:
1343	c
1344	REAL zgeom(klon,klev)
1345	REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
1346	REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
1347	REAL papmf(klon,klev),papmh(klon,klev+1)
1348	c
1349	c initialiser les variables de sortie (pour securite)
1350	c
1351	DO i = 1,klon
1352	pulow(i) = 0.0
1353	pvlow(i) = 0.0
1354	pustr(i) = 0.0
1355	pvstr(i) = 0.0
1356	ENDDO
1357	DO k = 1, klev
1358	DO i = 1, klon
1359	d_t(i,k) = 0.0
1360	d_u(i,k) = 0.0
1361	d_v(i,k) = 0.0
1362	pdudt(i,k)=0.0
1363	pdvdt(i,k)=0.0
1364	pdtdt(i,k)=0.0
1365	ENDDO
1366	ENDDO
1367	c
1368	c preparer les variables d'entree (attention: l'ordre des niveaux
1369	c verticaux augmente du haut vers le bas)
1370	c
1371	DO k = 1, klev
1372	DO i = 1, klon
1373	pt(i,k) = t(i,klev-k+1)
1374	pu(i,k) = u(i,klev-k+1)
1375	pv(i,k) = v(i,klev-k+1)
1376	papmf(i,k) = pplay(i,klev-k+1)
1377	ENDDO
1378	ENDDO
1379	DO k = 1, klev+1
1380	DO i = 1, klon
1381	papmh(i,k) = paprs(i,klev-k+2)
1382	ENDDO
1383	ENDDO
1384	DO i = 1, klon
1385	zgeom(i,klev) = RD * pt(i,klev)
1386	. * LOG(papmh(i,klev+1)/papmf(i,klev))
1387	ENDDO
1388	DO k = klev-1, 1, -1
1389	DO i = 1, klon
1390	zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
1391	. * LOG(papmf(i,k+1)/papmf(i,k))
1392	ENDDO
1393	ENDDO
1394	c
1395	c appeler la routine principale
1396	c
1397	CALL OROLIFT(klon,klev,ktest,
1398	. dtime,
1399	. papmh, zgeom,
1400	. pt, pu, pv,
1401	. plat,pmea, pstd, ppic,
1402	. pulow,pvlow,
1403	. pdudt,pdvdt,pdtdt)
1404	C
1405	DO k = 1, klev
1406	DO i = 1, klon
1407	d_u(i,klev+1-k) = dtime*pdudt(i,k)
1408	d_v(i,klev+1-k) = dtime*pdvdt(i,k)
1409	d_t(i,klev+1-k) = dtime*pdtdt(i,k)
1410	pustr(i) = pustr(i)
1411	. +RGpdudt(i,k)(papmh(i,k+1)-papmh(i,k))
1412	pvstr(i) = pvstr(i)
1413	. +RGpdvdt(i,k)(papmh(i,k+1)-papmh(i,k))
1414	ENDDO
1415	ENDDO
1416	c
1417	RETURN
1418	END
1419	SUBROUTINE OROLIFT( NLON,NLEV
1420	I , KTEST
1421	R , PTSPHY
1422	R , PAPHM1,PGEOM1,PTM1,PUM1,PVM1
1423	R , PLAT
1424	R , PMEA, PVAROR, ppic
1425	C OUTPUTS
1426	R , PULOW,PVLOW
1427	R , PVOM,PVOL,PTE )
1428
1429	C
1430	C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
1431	C
1432	C PURPOSE.
1433	C --------
1434	C
1435	C** INTERFACE.
1436	C ----------
1437	C CALLED FROM *lift_noro
1438	C ----------
1439	C
1440	C AUTHOR.
1441	C -------
1442	C F.LOTT LMD 22/11/95
1443	C
1444	USE dimphy
1445	implicit none
1446	C
1447	C
1448	cym#include "dimensions.h"
1449	cym#include "dimphy.h"
1450	#include "YOMCST.h"
1451	#include "YOEGWD.h"
1452	C-----------------------------------------------------------------------
1453	C
1454	C* 0.1 ARGUMENTS
1455	C ---------
1456	C
1457	C
1458	integer nlon, nlev
1459	REAL PTE(NLON,NLEV),
1460	* PVOL(NLON,NLEV),
1461	* PVOM(NLON,NLEV),
1462	* PULOW(NLON),
1463	* PVLOW(NLON)
1464	REAL PUM1(NLON,NLEV),
1465	* PVM1(NLON,NLEV),
1466	* PTM1(NLON,NLEV),
1467	* PLAT(NLON),PMEA(NLON),
1468	* PVAROR(NLON),
1469	* ppic(NLON),
1470	* PGEOM1(NLON,NLEV),
1471	* PAPHM1(NLON,NLEV+1)
1472	C
1473	INTEGER KTEST(NLON)
1474	real ptsphy
1475	C-----------------------------------------------------------------------
1476	C
1477	C* 0.2 LOCAL ARRAYS
1478	C ------------
1479	logical lifthigh, ll1
1480	cym integer klevm1, jl, ilevh, jk
1481	integer jl, ilevh, jk
1482	real zcons1, ztmst, zrtmst,zpi, zhgeo
1483	real zdelp, zslow, zsqua, zscav, zbet
1484	INTEGER
1485	* IKNUB(klon),
1486	* IKNUL(klon)
1487	LOGICAL LL1(KLON,KLEV+1)
1488	C
1489	REAL ZTAU(KLON,KLEV+1),
1490	* ZTAV(KLON,KLEV+1),
1491	* ZRHO(KLON,KLEV+1)
1492	REAL ZDUDT(KLON),
1493	* ZDVDT(KLON)
1494	REAL ZHCRIT(KLON,KLEV)
1495	C-----------------------------------------------------------------------
1496	C
1497	C* 1.1 INITIALIZATIONS
1498	C ---------------
1499
1500	LIFTHIGH=.FALSE.
1501
1502	IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)STOP
1503	ZCONS1=1./RD
1504	cym KLEVM1=KLEV-1
1505	ZTMST=PTSPHY
1506	ZRTMST=1./ZTMST
1507	ZPI=ACOS(-1.)
1508	C
1509	DO 1001 JL=kidia,kfdia
1510	ZRHO(JL,KLEV+1) =0.0
1511	PULOW(JL) =0.0
1512	PVLOW(JL) =0.0
1513	iknub(JL) =klev
1514	iknul(JL) =klev
1515	ilevh=klev/3
1516	ll1(jl,klev+1)=.false.
1517	DO 1000 JK=1,KLEV
1518	PVOM(JL,JK)=0.0
1519	PVOL(JL,JK)=0.0
1520	PTE (JL,JK)=0.0
1521	1000 CONTINUE
1522	1001 CONTINUE
1523
1524	C
1525	C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
1526	C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
1527	C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
1528	C
1529	C
1530	C
1531	DO 2006 JK=KLEV,1,-1
1532	DO 2007 JL=kidia,kfdia
1533	IF(KTEST(JL).EQ.1) THEN
1534	ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.)
1535	ZHGEO=PGEOM1(JL,JK)/RG
1536	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
1537	IF(ll1(JL,JK).XOR.ll1(JL,JK+1)) THEN
1538	iknub(JL)=JK
1539	ENDIF
1540	ENDIF
1541	2007 CONTINUE
1542	2006 CONTINUE
1543	C
1544	do 2010 jl=kidia,kfdia
1545	IF(KTEST(JL).EQ.1) THEN
1546	iknub(jl)=max(iknub(jl),klev/2)
1547	iknul(jl)=max(iknul(jl),2*klev/3)
1548	if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
1549	if(iknub(jl).eq.nktopg) iknul(jl)=klev
1550	if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
1551	ENDIF
1552	2010 continue
1553
1554	C do 2011 jl=kidia,kfdia
1555	C IF(KTEST(JL).EQ.1) THEN
1556	C print *,' iknul= ',iknul(jl),' iknub=',iknub(jl)
1557	C ENDIF
1558	C2011 continue
1559
1560	C PRINT *,' DANS OROLIFT: 2010'
1561
1562	DO 223 JK=KLEV,2,-1
1563	DO 222 JL=kidia,kfdia
1564	ZRHO(JL,JK)=2.PAPHM1(JL,JK)ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
1565	222 CONTINUE
1566	223 CONTINUE
1567	C PRINT *,' DANS OROLIFT: 223'
1568
1569	C********************************************************************
1570	C
1571	C* DEFINE LOW LEVEL FLOW
1572	C -------------------
1573	DO 2115 JK=klev,1,-1
1574	DO 2116 JL=kidia,kfdia
1575	IF(KTEST(JL).EQ.1) THEN
1576	if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
1577	pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1578	pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1579	zrho(JL,klev+1)=zrho(JL,klev+1)
1580	* +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1581	end if
1582	ENDIF
1583	2116 CONTINUE
1584	2115 CONTINUE
1585	DO 2110 JL=kidia,kfdia
1586	IF(KTEST(JL).EQ.1) THEN
1587	pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1588	pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1589	zrho(JL,klev+1)=zrho(JL,klev+1)
1590	* /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1591	ENDIF
1592	2110 CONTINUE
1593
1594
1595	200 CONTINUE
1596
1597	C***********************************************************
1598	C
1599	C* 3. COMPUTE MOUNTAIN LIFT
1600	C
1601	300 CONTINUE
1602	C
1603	DO 301 JL=kidia,kfdia
1604	IF(KTEST(JL).EQ.1) THEN
1605	ZTAU(JL,KLEV+1)= - GKLIFTZRHO(JL,KLEV+1)2.ROMEGA
1606	C * (2PVAROR(JL)+PMEA(JL))
1607	* 2PVAROR(JL)
1608	* SIN(ZPI/180.PLAT(JL))PVLOW(JL)
1609	ZTAV(JL,KLEV+1)= GKLIFTZRHO(JL,KLEV+1)2.ROMEGA
1610	C * (2PVAROR(JL)+PMEA(JL))
1611	* 2PVAROR(JL)
1612	* SIN(ZPI/180.PLAT(JL))PULOW(JL)
1613	ELSE
1614	ZTAU(JL,KLEV+1)=0.0
1615	ZTAV(JL,KLEV+1)=0.0
1616	ENDIF
1617	301 CONTINUE
1618
1619	C
1620	C* 4. COMPUTE LIFT PROFILE
1621	C* --------------------
1622	C
1623
1624	400 CONTINUE
1625
1626	DO 401 JK=1,KLEV
1627	DO 401 JL=kidia,kfdia
1628	IF(KTEST(JL).EQ.1) THEN
1629	ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
1630	ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
1631	ELSE
1632	ZTAU(JL,JK)=0.0
1633	ZTAV(JL,JK)=0.0
1634	ENDIF
1635	401 CONTINUE
1636	C
1637	C
1638	C* 5. COMPUTE TENDENCIES.
1639	C* -------------------
1640	IF(LIFTHIGH)THEN
1641	C
1642	500 CONTINUE
1643	C PRINT *,' DANS OROLIFT: 500'
1644	C
1645	C EXPLICIT SOLUTION AT ALL LEVELS
1646	C
1647	DO 524 JK=1,klev
1648	DO 523 JL=KIDIA,KFDIA
1649	IF(KTEST(JL).EQ.1) THEN
1650	ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
1651	ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP
1652	ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP
1653	ENDIF
1654	523 CONTINUE
1655	524 CONTINUE
1656	C
1657	C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
1658	C
1659	DO 530 JK=1,klev
1660	DO 530 JL=KIDIA,KFDIA
1661	IF(KTEST(JL).EQ.1) THEN
1662
1663	ZSLOW=SQRT(PULOW(JL)2+PVLOW(JL)2)
1664	ZSQUA=AMAX1(SQRT(PUM1(JL,JK)2+PVM1(JL,JK)2),GVSEC)
1665	ZSCAV=-ZDUDT(JL)PVM1(JL,JK)+ZDVDT(JL)PUM1(JL,JK)
1666	IF(ZSQUA.GT.GVSEC)THEN
1667	PVOM(JL,JK)=-ZSCAVPVM1(JL,JK)/ZSQUA*2
1668	PVOL(JL,JK)= ZSCAVPUM1(JL,JK)/ZSQUA*2
1669	ELSE
1670	PVOM(JL,JK)=0.0
1671	PVOL(JL,JK)=0.0
1672	ENDIF
1673	ZSQUA=SQRT(PUM1(JL,JK)2+PUM1(JL,JK)2)
1674	IF(ZSQUA.LT.ZSLOW)THEN
1675	PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK)
1676	PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK)
1677	ENDIF
1678
1679	ENDIF
1680	530 CONTINUE
1681	C
1682	C 6. LOW LEVEL LIFT, SEMI IMPLICIT:
1683	C ----------------------------------
1684
1685	ELSE
1686
1687	DO 601 JL=KIDIA,KFDIA
1688	IF(KTEST(JL).EQ.1) THEN
1689	DO JK=KLEV,IKNUB(JL),-1
1690	ZBET=GKLIFT2.ROMEGASIN(ZPI/180.PLAT(JL))ztmst
1691	* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL, JK))/
1692	* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV))
1693	ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2)
1694	ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2)
1695	PVOM(JL,JK)= ZBET*2ZDUDT(JL) - ZBET *ZDVDT(JL)
1696	PVOL(JL,JK)= ZBET ZDUDT(JL) + ZBET2ZDVDT(JL)
1697	ENDDO
1698	ENDIF
1699	601 CONTINUE
1700
1701	ENDIF
1702
1703	RETURN
1704	END
1705	SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay)
1706	C
1707	C**** SUGWD INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
1708	C
1709	C PURPOSE.
1710	C --------
1711	C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
1712	C GRAVITY WAVE DRAG PARAMETRIZATION.
1713	C
1714	C** INTERFACE.
1715	C ----------
1716	C CALL SUGWD FROM SUPHEC
1717	C ----- ------
1718	C
1719	C EXPLICIT ARGUMENTS :
1720	C --------------------
1721	C PSIG : VERTICAL COORDINATE TABLE
1722	C NLEV : NUMBER OF MODEL LEVELS
1723	C
1724	C IMPLICIT ARGUMENTS :
1725	C --------------------
1726	C COMMON YOEGWD
1727	C
1728	C METHOD.
1729	C -------
1730	C SEE DOCUMENTATION
1731	C
1732	C EXTERNALS.
1733	C ----------
1734	C NONE
1735	C
1736	C REFERENCE.
1737	C ----------
1738	C ECMWF Research Department documentation of the IFS
1739	C
1740	C AUTHOR.
1741	C -------
1742	C MARTIN MILLER ECMWF
1743	C
1744	C MODIFICATIONS.
1745	C --------------
1746	C ORIGINAL : 90-01-01
1747	C ------------------------------------------------------------------
1748	implicit none
1749	C
1750	C -----------------------------------------------------------------
1751	#include "YOEGWD.h"
1752	C ----------------------------------------------------------------
1753	C
1754	integer nlon,nlev, jk
1755	REAL paprs(nlon,nlev+1)
1756	REAL pplay(nlon,nlev)
1757	real zpr,zstra,zsigt,zpm1r
1758	C
1759	C* 1. SET THE VALUES OF THE PARAMETERS
1760	C --------------------------------
1761	C
1762	100 CONTINUE
1763	C
1764	PRINT *,' DANS SUGWD NLEV=',NLEV
1765	GHMAX=10000.
1766	C
1767	ZPR=100000.
1768	ZSTRA=0.1
1769	ZSIGT=0.94
1770	cold ZPR=80000.
1771	cold ZSIGT=0.85
1772	C
1773	DO 110 JK=1,NLEV
1774	ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1)
1775	IF(ZPM1R.GE.ZSIGT)THEN
1776	nktopg=JK
1777	ENDIF
1778	ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1)
1779	IF(ZPM1R.GE.ZSTRA)THEN
1780	NSTRA=JK
1781	ENDIF
1782	110 CONTINUE
1783	c
1784	c inversion car dans orodrag on compte les niveaux a l'envers
1785	nktopg=nlev-nktopg+1
1786	nstra=nlev-nstra
1787	print *,' DANS SUGWD nktopg=', nktopg
1788	print *,' DANS SUGWD nstra=', nstra
1789	C
1790	GSIGCR=0.80
1791	C
1792	GKDRAG=0.2
1793	GRAHILO=1.
1794	GRCRIT=0.01
1795	GFRCRIT=1.0
1796	GKWAKE=0.50
1797	C
1798	GKLIFT=0.50
1799	GVCRIT =0.0
1800	C
1801	C
1802	C ----------------------------------------------------------------
1803	C
1804	C* 2. SET VALUES OF SECURITY PARAMETERS
1805	C ---------------------------------
1806	C
1807	200 CONTINUE
1808	C
1809	GVSEC=0.10
1810	GSSEC=1.E-12
1811	C
1812	GTSEC=1.E-07
1813	C
1814	C ----------------------------------------------------------------
1815	C
1816	RETURN
1817	END

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

Download in other formats:

Original Format