Context Navigation

orografi.F @ 1990

Last change on this file since 1990 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 48.5 KB

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

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

Download in other formats:

Original Format