Context Navigation

orosetup.F @ 3303

Last change on this file since 3303 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

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

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

Download in other formats:

Original Format