Context Navigation

orosetup.F @ 1243

Last change on this file since 1243 was 1056, checked in by slebonnois, 11 years ago
SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 15.0 KB

Line
1	SUBROUTINE orosetup
2	* ( nlon , nlev , ktest
3	* , kkcrit, kkcrith, kcrit, ksect , kkhlim
4	* , kkenvh, kknu , kknu2
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 "dimensions.h"
101	#include "paramet.h"
102
103	#include "YOMCST.h"
104	#include "YOEGWD.h"
105
106	c-----------------------------------------------------------------------
107	c
108	c* 0.1 arguments
109	c ---------
110	c
111	integer nlon,nlev
112	integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
113	* kkhlim(nlon),ktest(nlon),kkenvh(nlon)
114
115	c
116	real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
117	* pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
118	* prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
119	* ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
120	* pzdep(nlon,klev)
121	real pulow(nlon),pvlow(nlon),ptheta(nlon),pgam(nlon),pnu(nlon),
122	* pd1(nlon),pd2(nlon),pdmod(nlon)
123	real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
124	c
125	c-----------------------------------------------------------------------
126	c
127	c* 0.2 local arrays
128	c ------------
129	c
130	c
131	integer ilevh ,jl,jk,iii
132	real zcons1,zhgeo,zu,zphi
133	real zvt1,zvt2,zdwind,zwind,zdelp
134	real zstabm,zstabp,zrhom,zrhop
135	logical lo
136	logical ll1(klon,klev+1)
137	integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
138	* kentp(klon),ncount(klon)
139	c
140	real zhcrit(klon,klev),zvpf(klon,klev),
141	* zdp(klon,klev)
142	real znorm(klon),zb(klon),zc(klon),
143	* zulow(klon),zvlow(klon),znup(klon),znum(klon)
144
145	c ------------------------------------------------------------------
146	c
147	c* 1. initialization
148	c --------------
149	c
150	c PRINT *,' in orosetup'
151	100 continue
152	c
153	c ------------------------------------------------------------------
154	c
155	c* 1.1 computational constants
156	c -----------------------
157	c
158	110 continue
159	c
160	ilevh =klev/3
161	c
162	zcons1=1./rd
163	c
164	c ------------------------------------------------------------------
165	c
166	c* 2.
167	c --------------
168	c
169	200 continue
170	c
171	c ------------------------------------------------------------------
172	c
173	c* 2.1 define low level wind, project winds in plane of
174	c* low level wind, determine sector in which to take
175	c* the variance and set indicator for critical levels.
176	c
177	c
178	c
179	do 2001 jl=kidia,kfdia
180	kknu(jl) =klev
181	kknu2(jl) =klev
182	kknub(jl) =klev
183	kknul(jl) =klev
184	pgam(jl) =max(pgam(jl),gtsec)
185	ll1(jl,klev+1)=.false.
186	2001 continue
187	c
188	c Ajouter une initialisation (L. Li, le 23fev99):
189	c
190	do jk=klev,ilevh,-1
191	do jl=kidia,kfdia
192	ll1(jl,jk)= .false.
193	ENDDO
194	ENDDO
195	c
196	c* define top of low level flow
197	c ----------------------------
198	do 2002 jk=klev,ilevh,-1
199	do 2003 jl=kidia,kfdia
200	if(ktest(jl).eq.1) then
201	lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
202	if(lo) then
203	kkcrit(jl)=jk
204	endif
205	zhcrit(jl,jk)=ppic(jl)-pval(jl)
206	zhgeo=pgeom1(jl,jk)/rg
207	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
208	C if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
209	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
210	kknu(jl)=jk
211	endif
212	if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
213	endif
214	2003 continue
215	2002 continue
216	do 2004 jk=klev,ilevh,-1
217	do 2005 jl=kidia,kfdia
218	if(ktest(jl).eq.1) then
219	zhcrit(jl,jk)=ppic(jl)-pmea(jl)
220	zhgeo=pgeom1(jl,jk)/rg
221	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
222	if(ll1(jl,jk) .neqv. ll1(jl,jk+1)) then
223	kknu2(jl)=jk
224	endif
225	if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
226	endif
227	2005 continue
228	2004 continue
229	do 2006 jk=klev,ilevh,-1
230	do 2007 jl=kidia,kfdia
231	if(ktest(jl).eq.1) then
232	zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
233	zhgeo=pgeom1(jl,jk)/rg
234	ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
235	c if(ll1(jl,jk).xor.ll1(jl,jk+1)) then
236	if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
237	kknub(jl)=jk
238	endif
239	if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
240	endif
241	2007 continue
242	2006 continue
243	c
244	do 2010 jl=kidia,kfdia
245	if(ktest(jl).eq.1) then
246	kknu(jl)=min(kknu(jl),nktopg)
247	kknu2(jl)=min(kknu2(jl),nktopg)
248	kknub(jl)=min(kknub(jl),nktopg)
249	kknul(jl)=klev
250	endif
251	2010 continue
252	c
253	210 continue
254	c
255	cc* initialize various arrays
256	c
257	do 2107 jl=kidia,kfdia
258	prho(jl,klev+1) =0.0
259	cym correction en attendant mieux
260	prho(jl,1) =0.0
261	pstab(jl,klev+1) =0.0
262	pstab(jl,1) =0.0
263	pri(jl,klev+1) =9999.0
264	ppsi(jl,klev+1) =0.0
265	pri(jl,1) =0.0
266	pvph(jl,1) =0.0
267	pvph(jl,klev+1) =0.0
268	cym correction en attendant mieux
269	cym pvph(jl,klev) =0.0
270	pulow(jl) =0.0
271	pvlow(jl) =0.0
272	zulow(jl) =0.0
273	zvlow(jl) =0.0
274	kkcrith(jl) =klev
275	kkenvh(jl) =klev
276	kentp(jl) =klev
277	kcrit(jl) =1
278	ncount(jl) =0
279	ll1(jl,klev+1) =.false.
280	2107 continue
281	c
282	c* define flow density and stratification (rho and N2)
283	c at semi layers.
284	c -------------------------------------------------------
285	c
286	do 223 jk=klev,2,-1
287	do 222 jl=kidia,kfdia
288	if(ktest(jl).eq.1) then
289	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
290	prho(jl,jk)=2.paphm1(jl,jk)zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
291	endif
292	222 continue
293	223 continue
294	c print*,"altitude(m)=",pgeom1(kfdia/2,:)
295	c print*,"pression(Pa)=",papm1(kfdia/2,:)
296	c
297	c********************************************************************
298	c
299	c* define Low level flow (between ground and peacks-valleys)
300	c ---------------------------------------------------------
301	do 2115 jk=klev,ilevh,-1
302	do 2116 jl=kidia,kfdia
303	if(ktest(jl).eq.1) then
304	if(jk.ge.kknu2(jl).and.jk.le.kknul(jl)) then
305	pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
306	pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
307	pstab(jl,klev+1)=pstab(jl,klev+1)
308	c +pstab(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
309	prho(jl,klev+1)=prho(jl,klev+1)
310	c +prho(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
311	end if
312	endif
313	2116 continue
314	2115 continue
315	do 2110 jl=kidia,kfdia
316	if(ktest(jl).eq.1) then
317	pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
318	pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
319	znorm(jl)=max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
320	pvph(jl,klev+1)=znorm(jl)
321	pstab(jl,klev+1)=pstab(jl,klev+1)
322	c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
323	prho(jl,klev+1)=prho(jl,klev+1)
324	c /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
325	endif
326	2110 continue
327
328	c
329	c******* setup orography orientation relative to the low level
330	C wind and define parameters of the Anisotropic wave stress.
331	c
332	do 2112 jl=kidia,kfdia
333	if(ktest(jl).eq.1) then
334	lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
335	if(lo) then
336	zu=pulow(jl)+2.*gvsec
337	else
338	zu=pulow(jl)
339	endif
340	zphi=atan(pvlow(jl)/zu)
341	ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
342	zb(jl)=1.-0.18pgam(jl)-0.04pgam(jl)**2
343	zc(jl)=0.48pgam(jl)+0.3pgam(jl)**2
344	pd1(jl)=zb(jl)-(zb(jl)-zc(jl))(sin(ppsi(jl,klev+1))*2)
345	pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))
346	* *cos(ppsi(jl,klev+1))
347	pdmod(jl)=sqrt(pd1(jl)2+pd2(jl)2)
348	endif
349	2112 continue
350	c
351	c ********** projet flow in plane of lowlevel stress ***********
352	C ********** Find critical levels... ***********
353	c
354	do 213 jk=1,klev
355	do 212 jl=kidia,kfdia
356	if(ktest(jl).eq.1) then
357	zvt1 =pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk)
358	zvt2 =-pvlow(jl)pum1(jl,jk)+pulow(jl)pvm1(jl,jk)
359	zvpf(jl,jk)=(zvt1pd1(jl)+zvt2pd2(jl))/(znorm(jl)*pdmod(jl))
360	endif
361	ptau(jl,jk) =0.0
362	pzdep(jl,jk) =0.0
363	ppsi(jl,jk) =0.0
364	ll1(jl,jk) =.false.
365	212 continue
366	213 continue
367	do 215 jk=2,klev
368	do 214 jl=kidia,kfdia
369	if(ktest(jl).eq.1) then
370	zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
371	pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
372	* (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
373	* /zdp(jl,jk)
374	if(pvph(jl,jk).lt.gvsec) then
375	pvph(jl,jk)=gvsec
376	kcrit(jl)=jk
377	endif
378	endif
379	214 continue
380	215 continue
381	c
382	c* 2.3 mean flow richardson number.
383	c
384	230 continue
385	c
386	do 232 jk=2,klev
387	do 231 jl=kidia,kfdia
388	if(ktest(jl).eq.1) then
389	zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
390	pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
391	* /(rgprho(jl,jk)zdwind))**2
392	pri(jl,jk)=max(pri(jl,jk),grcrit)
393	endif
394	231 continue
395	232 continue
396
397	c
398	c
399	c* define top of 'envelope' layer
400	c ----------------------------
401
402	do 233 jl=kidia,kfdia
403	pnu (jl)=0.0
404	znum(jl)=0.0
405	233 continue
406
407	do 234 jk=2,klev-1
408	do 234 jl=kidia,kfdia
409
410	if(ktest(jl).eq.1) then
411
412	if (jk.ge.kknu2(jl)) then
413
414	znum(jl)=pnu(jl)
415	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
416	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
417	zwind=max(sqrt(zwind**2),gvsec)
418	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
419	zstabm=sqrt(max(pstab(jl,jk ),gssec))
420	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
421	zrhom=prho(jl,jk )
422	zrhop=prho(jl,jk+1)
423	pnu(jl) = pnu(jl) + (zdelp/rg)*
424	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
425	if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
426	* .and.(kkenvh(jl).eq.klev))
427	* kkenvh(jl)=jk
428
429	endif
430
431	endif
432
433	234 continue
434
435	c calculation of a dynamical mixing height for when the waves
436	C BREAK AT LOW LEVEL: The drag will be repartited over
437	C a depths that depends on waves vertical wavelength,
438	C not just between two adjacent model layers.
439	c of gravity waves:
440
441	do 235 jl=kidia,kfdia
442	znup(jl)=0.0
443	znum(jl)=0.0
444	235 continue
445
446	do 236 jk=klev-1,2,-1
447	do 236 jl=kidia,kfdia
448
449	if(ktest(jl).eq.1) then
450
451	znum(jl)=znup(jl)
452	zwind=(pulow(jl)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
453	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
454	zwind=max(sqrt(zwind**2),gvsec)
455	zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
456	zstabm=sqrt(max(pstab(jl,jk ),gssec))
457	zstabp=sqrt(max(pstab(jl,jk+1),gssec))
458	zrhom=prho(jl,jk )
459	zrhop=prho(jl,jk+1)
460	znup(jl) = znup(jl) + (zdelp/rg)*
461	* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
462	if((znum(jl).le.rpi/4.).and.(znup(jl).gt.rpi/4.)
463	* .and.(kkcrith(jl).eq.klev))
464	* kkcrith(jl)=jk
465
466	endif
467
468	236 continue
469
470	do 237 jl=kidia,kfdia
471	if(ktest(jl).eq.1) then
472	kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
473	kkcrith(jl)=max0(kkcrith(jl),kknu(jl))
474	if(kcrit(jl).ge.kkcrith(jl))kcrit(jl)=1
475	endif
476	237 continue
477	c
478	c directional info for flow blocking *************************
479	c
480	do 251 jk=1,klev
481	do 252 jl=kidia,kfdia
482	if(ktest(jl).eq.1) then
483	lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
484	if(lo) then
485	zu=pum1(jl,jk)+2.*gvsec
486	else
487	zu=pum1(jl,jk)
488	endif
489	zphi=atan(pvm1(jl,jk)/zu)
490	ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
491	endif
492	252 continue
493	251 continue
494
495	c forms the vertical 'leakiness' **************************
496
497	do 254 jk=ilevh,klev
498	do 253 jl=kidia,kfdia
499	if(ktest(jl).eq.1) then
500	pzdep(jl,jk)=0
501	if(jk.ge.kkenvh(jl).and.kkenvh(jl).ne.klev) then
502	pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl) )-pgeom1(jl, jk))/
503	* (pgeom1(jl,kkenvh(jl) )-pgeom1(jl,klev))
504	end if
505	endif
506	253 continue
507	254 continue
508
509	return
510	end
511
512

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

Download in other formats:

Original Format