Context Navigation

orosetup.F @ 3884

Last change on this file since 3884 was 3884, checked in by ikovalenko, 4 months ago

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

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

Download in other formats:

Original Format