Context Navigation

orosetup.F @ 3461

Last change on this file since 3461 was 2047, checked in by slebonnois, 6 years ago
SL: VENUS, ajout des modifs apportees par Thomas Navarro pour la parametrisation des ondes de gravite orographiques
File size: 16.2 KB

Line
1	SUBROUTINE orosetup
2	* ( nlon , nlev , ktest
3	* , kkcrit, kkcrith, kcrit, ksect , kkhlim
4	* , kkenvh, kknu , kknu2, kkbreak
5	* , paphm1, papm1 , pum1 , pvm1, ptm1, pgeom1, pstab, pstd
6	* , prho , pri , ptau, pvph, ppsi, pzdep
7	* , pulow , pvlow
8	* , ptheta, pgam, pmea, ppic, pval
9	* , pnu , pd1 , pd2 ,pdmod )
10	C
11	c**** gwsetup
12	c
13	c purpose.
14	c --------
15	c SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME:
16	C DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND
17	C STRATIFICATION.....
18	c
19	c** interface.
20	c ----------
21	c from orodrag
22	c
23	c explicit arguments :
24	c --------------------
25	c ==== inputs ===
26	c
27	c nlon----input-I-Total number of horizontal points that get into physics
28	c nlev----input-I-Number of vertical levels
29	c ktest--input-I: Flags to indicate active points
30	c
31	c ptsphy--input-R-Time-step (s)
32	c paphm1--input-R: pressure at model 1/2 layer
33	c papm1---input-R: pressure at model layer
34	c pgeom1--input-R: Altitude of layer above ground
35	c VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
36	c pstab-----R-: Brunt-Vaisala freq.^2 at 1/2 layers (input except klev+1)
37	c ptm1, pum1, pvm1--R-: t, u and v
38	c pmea----input-R-Mean Orography (m)
39	C pstd----input-R-SSO standard deviation (m)
40	c psig----input-R-SSO slope
41	c pgam----input-R-SSO Anisotropy
42	c pthe----input-R-SSO Angle
43	c ppic----input-R-SSO Peacks elevation (m)
44	c pval----input-R-SSO Valleys elevation (m)
45
46	c ==== outputs ===
47	c pulow, pvlow -output-R: Low-level wind
48	c kkcrit----I-: Security value for top of low level flow
49	c kcrit-----I-: Critical level
50	c ksect-----I-: Not used
51	c kkhlim----I-: Not used
52	c kkenvh----I-: Top of blocked flow layer
53	c kknu------I-: Layer that sees mountain peacks
54	c kknu2-----I-: Layer that sees mountain peacks above mountain mean
55	c kknub-----I-: Layer that sees mountain mean above valleys
56	c prho------R-: Density at 1/2 layers
57	c pri-------R-: Background Richardson Number, Wind shear measured along GW stress
58	c pvph------R-: Wind in plan of GW stress, Half levels.
59	c ppsi------R-: Angle between low level wind and SS0 main axis.
60	c pd1-------R-\| Compared the ratio of the stress
61	c pd2-------R-\| that is along the wind to that Normal to it.
62	c pdi define the plane of low level stress
63	c compared to the low level wind.
64	c see p. 108 Lott & Miller (1997).
65	c pdmod-----R-: Norme of pdi
66
67	c === local arrays ===
68	c
69	c zvpf------R-: Wind projected in the plan of the low-level stress.
70
71	c ==== outputs ===
72	c
73	c implicit arguments : none
74	c --------------------
75	c
76	c method.
77	c -------
78	c
79	c
80	c externals.
81	c ----------
82	c
83	c
84	c reference.
85	c ----------
86	c
87	c see ecmwf research department documentation of the "i.f.s."
88	c
89	c author.
90	c -------
91	c
92	c modifications.
93	c --------------
94	c f.lott for the new-gwdrag scheme november 1993
95	c
96	c-----------------------------------------------------------------------
97	use dimphy
98	implicit none
99
100	#include "YOMCST.h"
101	#include "YOEGWD.h"
102
103	c-----------------------------------------------------------------------
104	c
105	c* 0.1 arguments
106	c ---------
107	c
108	integer nlon,nlev
109	integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
110	* kkhlim(nlon),ktest(nlon),kkenvh(nlon)
111	integer kkbreak(nlon)
112	c
113	real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
114	* pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
115	* prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
116	* ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
117	* pzdep(nlon,klev)
118	real pulow(nlon),pvlow(nlon),ptheta(nlon),pgam(nlon),pnu(nlon),
119	* pd1(nlon),pd2(nlon),pdmod(nlon)
120	real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
121	c
122	c-----------------------------------------------------------------------
123	c
124	c* 0.2 local arrays
125	c ------------
126	c
127	c
128	integer ilevh ,jl,jk,iii
129	real zcons1,zhgeo,zu,zphi
130	real zvt1,zvt2,zdwind,zwind,zdelp
131	real zstabm,zstabp,zrhom,zrhop
132	logical lo
133	logical ll1(klon,klev+1)
134	integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
135	* kentp(klon),ncount(klon)
136	c
137	real zhcrit(klon,klev),zvpf(klon,klev),
138	* zdp(klon,klev)
139	real znorm(klon),zb(klon),zc(klon),
140	* zulow(klon),zvlow(klon),znup(klon),znum(klon)
141
142	c ------------------------------------------------------------------
143	c
144	c* 1. initialization
145	c --------------
146	c
147	c PRINT *,' in orosetup'
148	100 continue
149	c
150	c ------------------------------------------------------------------
151	c
152	c* 1.1 computational constants
153	c -----------------------
154	c
155	110 continue
156	c
157	ilevh =klev/3
158	c
159	zcons1=1./rd
160	c
161	c ------------------------------------------------------------------
162	c
163	c* 2.
164	c --------------
165	c
166	200 continue
167	c
168	c ------------------------------------------------------------------
169	c
170	c* 2.1 define low level wind, project winds in plane of
171	c* low level wind, determine sector in which to take
172	c* the variance and set indicator for critical levels.
173	c
174	c
175	c
176	do 2001 jl=kidia,kfdia
177	kknu(jl) =klev
178	kknu2(jl) =klev
179	kknub(jl) =klev
180	kknul(jl) =klev
181	kkbreak(jl) =klev + 1
182	pgam(jl) =max(pgam(jl),gtsec)
183	ll1(jl,klev+1)=.false.
184	2001 continue
185	c
186	c Ajouter une initialisation (L. Li, le 23fev99):
187	c
188	do jk=klev,ilevh,-1
189	do jl=kidia,kfdia
190	ll1(jl,jk)= .false.
191	ENDDO
192	ENDDO
193
194	c VENUS: define break for subcritical stress
195	c ----------------------------
196	do jk=klev,ilevh,-1
197	do jl=kidia,kfdia
198	if(ktest(jl).eq.1) then
199	!zhgeo=pgeom1(jl,jk)/rg
200	!!zhcrit(jl,jk)=ppic(jl)
201	!zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
202	!ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
203	!if(ll1(jl,jk) .neqv. ll1(jl,jk+1)) then
204	! kkbreak(jl)=jk
205	!endif
206
207	!if (paphm1(jl,jk) .ge. 7.e6) kkbreak(jl)=jk
208	!kkbreak(jl) = klev - 2 ! gwd1103
209	!kkbreak(jl) = klev - 4 ! gwd1104
210	!kkbreak(jl) = klev - 3 ! gwd1105
211
212	endif
213	enddo
214	enddo
215
216	c
217	c* define top of low level flow
218	c ----------------------------
219	do 2002 jk=klev,ilevh,-1
220	do 2003 jl=kidia,kfdia
221	if(ktest(jl).eq.1) then
222	lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
223	if(lo) then
224	kkcrit(jl)=jk
225	endif
226	zhcrit(jl,jk)=ppic(jl)-pval(jl)
227	zhgeo=pgeom1(jl,jk)/rg
228	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
229	C if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
230	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
231	kknu(jl)=jk
232	endif
233	if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
234	endif
235	2003 continue
236	2002 continue
237	do 2004 jk=klev,ilevh,-1
238	do 2005 jl=kidia,kfdia
239	if(ktest(jl).eq.1) then
240	! zhcrit(jl,jk)=ppic(jl)-pmea(jl)
241	zhgeo=pgeom1(jl,jk)/rg
242	! ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
243	ll1(jl,jk)=(zhgeo.gt.0.5*pstd(jl)) ! TN : do not consider outlier peaks
244	! ll1(jl,jk)=(zhgeo.gt.pstd(jl)) ! TN : do not consider outlier peaks
245	! ll1(jl,jk)=(zhgeo.gt.2*pstd(jl)) ! TN : do not consider outlier peaks
246	if(ll1(jl,jk) .neqv. ll1(jl,jk+1)) then
247	kknu2(jl)=jk
248	endif
249	if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
250	endif
251	2005 continue
252	2004 continue
253
254	! do 2104 jk=klev,ilevh,-1
255	! do 2105 jl=kidia,kfdia
256	! if(ktest(jl).eq.1) then
257	! zhgeo=pgeom1(jl,jk)/rg
258	! ll1(jl,jk)=(zhgeo.gt.pstd(jl)) ! TN : do not consider outlier peaks
259	! if(ll1(jl,jk) .neqv. ll1(jl,jk+1)) then
260	! kknul(jl)=jk
261	! endif
262	! if(.not.ll1(jl,ilevh))kknul(jl)=ilevh
263	! endif
264	! 2105 continue
265	! 2104 continue
266
267
268	do 2006 jk=klev,ilevh,-1
269	do 2007 jl=kidia,kfdia
270	if(ktest(jl).eq.1) then
271	zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
272	zhgeo=pgeom1(jl,jk)/rg
273	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
274	c if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
275	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
276	kknub(jl)=jk
277	endif
278	if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
279	endif
280	2007 continue
281	2006 continue
282	c
283	do 2010 jl=kidia,kfdia
284	if(ktest(jl).eq.1) then
285	kknu(jl)=min(kknu(jl),nktopg)
286	kknu2(jl)=min(kknu2(jl),nktopg)
287	kknub(jl)=min(kknub(jl),nktopg)
288	! kknul(jl)=klev
289	endif
290	2010 continue
291	c
292	210 continue
293	c
294	cc* initialize various arrays
295	c
296	do 2107 jl=kidia,kfdia
297	prho(jl,klev+1) =0.0
298	cym correction en attendant mieux
299	prho(jl,1) =0.0
300	pstab(jl,klev+1) =0.0
301	pstab(jl,1) =0.0
302	pri(jl,klev+1) =9999.0
303	ppsi(jl,klev+1) =0.0
304	pri(jl,1) =0.0
305	pvph(jl,1) =0.0
306	pvph(jl,klev+1) =0.0
307	cym correction en attendant mieux
308	cym pvph(jl,klev) =0.0
309	pulow(jl) =0.0
310	pvlow(jl) =0.0
311	zulow(jl) =0.0
312	zvlow(jl) =0.0
313	kkcrith(jl) =klev
314	kkenvh(jl) =klev
315	kentp(jl) =klev
316	kcrit(jl) =1
317	ncount(jl) =0
318	ll1(jl,klev+1) =.false.
319	2107 continue
320	c
321	c* define flow density and stratification (rho and N2)
322	c at semi layers.
323	c -------------------------------------------------------
324	c
325	do 223 jk=klev,2,-1
326	do 222 jl=kidia,kfdia
327	if(ktest(jl).eq.1) then
328	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
329	prho(jl,jk)=2.paphm1(jl,jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
330	endif
331	222 continue
332	223 continue
333	c print*,"altitude(m)=",pgeom1(kfdia/2,:)
334	c print*,"pression(Pa)=",papm1(kfdia/2,:)
335	c
336	c********************************************************************
337	c
338	c* define Low level flow (between ground and peacks-valleys)
339	c ---------------------------------------------------------
340	do 2115 jk=klev,ilevh,-1
341	do 2116 jl=kidia,kfdia
342	if(ktest(jl).eq.1) then
343	if(jk.ge.kknu2(jl).and.jk.le.kknul(jl)) then
344	pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
345	pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
346	pstab(jl,klev+1)=pstab(jl,klev+1)
347	c +pstab(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
348	prho(jl,klev+1)=prho(jl,klev+1)
349	c +prho(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
350	end if
351	endif
352	2116 continue
353	2115 continue
354	do 2110 jl=kidia,kfdia
355	if(ktest(jl).eq.1) then
356	pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
357	pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
358	znorm(jl)=max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
359	pvph(jl,klev+1)=znorm(jl)
360	pstab(jl,klev+1)=pstab(jl,klev+1)
361	c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
362	prho(jl,klev+1)=prho(jl,klev+1)
363	c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
364	endif
365	2110 continue
366
367	c
368	c******* setup orography orientation relative to the low level
369	C wind and define parameters of the Anisotropic wave stress.
370	c
371	do 2112 jl=kidia,kfdia
372	if(ktest(jl).eq.1) then
373	lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
374	if(lo) then
375	zu=pulow(jl)+2.*gvsec
376	else
377	zu=pulow(jl)
378	endif
379	zphi=atan(pvlow(jl)/zu)
380	ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
381	zb(jl)=1.-0.18pgam(jl)-0.04pgam(jl)**2
382	zc(jl)=0.48pgam(jl)+0.3pgam(jl)**2
383	pd1(jl)=zb(jl)-(zb(jl)-zc(jl))(sin(ppsi(jl,klev+1))*2)
384	pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))
385	* *cos(ppsi(jl,klev+1))
386	pdmod(jl)=sqrt(pd1(jl)2+pd2(jl)2)
387	endif
388	2112 continue
389	c
390	c ********** projet flow in plane of lowlevel stress ***********
391	C ********** Find critical levels... ***********
392	c
393	do 213 jk=1,klev
394	do 212 jl=kidia,kfdia
395	if(ktest(jl).eq.1) then
396	zvt1 =pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk)
397	zvt2 =-pvlow(jl)pum1(jl,jk)+pulow(jl)pvm1(jl,jk)
398	zvpf(jl,jk)=(zvt1pd1(jl)+zvt2pd2(jl))/(znorm(jl)*pdmod(jl))
399	endif
400	ptau(jl,jk) =0.0
401	pzdep(jl,jk) =0.0
402	ppsi(jl,jk) =0.0
403	ll1(jl,jk) =.false.
404	212 continue
405	213 continue
406	do 215 jk=2,klev
407	do 214 jl=kidia,kfdia
408	if(ktest(jl).eq.1) then
409	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
410	pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
411	* (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
412	* /zdp(jl,jk)
413	if(pvph(jl,jk).lt.gvsec) then
414	pvph(jl,jk)=gvsec
415	kcrit(jl)=jk
416	endif
417	endif
418	214 continue
419	215 continue
420	c
421	c* 2.3 mean flow richardson number.
422	c
423	230 continue
424	c
425	do 232 jk=2,klev
426	do 231 jl=kidia,kfdia
427	if(ktest(jl).eq.1) then
428	zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
429	pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
430	* /(rgprho(jl,jk)zdwind))**2
431	pri(jl,jk)=max(pri(jl,jk),grcrit)
432	endif
433	231 continue
434	232 continue
435
436	c
437	c
438	c* define top of 'envelope' layer
439	c ----------------------------
440
441	do 233 jl=kidia,kfdia
442	pnu (jl)=0.0
443	znum(jl)=0.0
444	233 continue
445
446	do 234 jk=2,klev-1
447	do 234 jl=kidia,kfdia
448
449	if(ktest(jl).eq.1) then
450
451	if (jk.ge.kknu2(jl)) then
452
453	znum(jl)=pnu(jl)
454	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
455	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
456	zwind=max(sqrt(zwind**2),gvsec)
457	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
458	zstabm=sqrt(max(pstab(jl,jk ),gssec))
459	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
460	zrhom=prho(jl,jk )
461	zrhop=prho(jl,jk+1)
462	pnu(jl) = pnu(jl) + (zdelp/rg)*
463	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
464	if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
465	* .and.(kkenvh(jl).eq.klev))
466	* kkenvh(jl)=jk
467
468	endif
469
470	endif
471
472	234 continue
473
474	c calculation of a dynamical mixing height for when the waves
475	C BREAK AT LOW LEVEL: The drag will be repartited over
476	C a depths that depends on waves vertical wavelength,
477	C not just between two adjacent model layers.
478	c of gravity waves:
479
480	do 235 jl=kidia,kfdia
481	znup(jl)=0.0
482	znum(jl)=0.0
483	235 continue
484
485	do 236 jk=klev-1,2,-1
486	do 236 jl=kidia,kfdia
487
488	if(ktest(jl).eq.1) then
489
490	znum(jl)=znup(jl)
491	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
492	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
493	zwind=max(sqrt(zwind**2),gvsec)
494	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
495	zstabm=sqrt(max(pstab(jl,jk ),gssec))
496	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
497	zrhom=prho(jl,jk )
498	zrhop=prho(jl,jk+1)
499	znup(jl) = znup(jl) + (zdelp/rg)*
500	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
501	if((znum(jl).le.rpi/4.).and.(znup(jl).gt.rpi/4.)
502	* .and.(kkcrith(jl).eq.klev))
503	* kkcrith(jl)=jk
504
505	endif
506
507	236 continue
508
509	do 237 jl=kidia,kfdia
510	if(ktest(jl).eq.1) then
511	kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
512	kkcrith(jl)=max0(kkcrith(jl),kknu(jl))
513	if(kcrit(jl).ge.kkcrith(jl))kcrit(jl)=1
514	endif
515	237 continue
516	c
517	c directional info for flow blocking *************************
518	c
519	do 251 jk=1,klev
520	do 252 jl=kidia,kfdia
521	if(ktest(jl).eq.1) then
522	lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
523	if(lo) then
524	zu=pum1(jl,jk)+2.*gvsec
525	else
526	zu=pum1(jl,jk)
527	endif
528	zphi=atan(pvm1(jl,jk)/zu)
529	ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
530	endif
531	252 continue
532	251 continue
533
534	c forms the vertical 'leakiness' **************************
535
536	do 254 jk=ilevh,klev
537	do 253 jl=kidia,kfdia
538	if(ktest(jl).eq.1) then
539	pzdep(jl,jk)=0
540	if(jk.ge.kkenvh(jl).and.kkenvh(jl).ne.klev) then
541	pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl) )-pgeom1(jl, jk))/
542	* (pgeom1(jl,kkenvh(jl) )-pgeom1(jl,klev))
543	end if
544	endif
545	253 continue
546	254 continue
547
548	return
549	end
550
551

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

Download in other formats:

Original Format