Context Navigation

orosetup.F @ 2156

Last change on this file since 2156 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

File size: 13.6 KB

Line
1	SUBROUTINE OROSETUP
2	* ( klon , klev , KTEST
3	* , KKCRIT, KKCRITH, KCRIT, KSECT , KKHLIM
4	* , kkenvh, kknu , kknu2
5	* , PAPHM1, PAPM1 , PUM1 , PVM1 , PTM1 , PGEOM1, pvaror
6	* , PRHO , PRI , PSTAB , PTAU , PVPH ,ppsi, pzdep
7	* , PULOW , PVLOW
8	* , Ptheta, pgamma, pnu , pd1 , pd2 ,pdmod )
9	C
10	C**** GWSETUP
11	C
12	C PURPOSE.
13	C --------
14	C
15	C** INTERFACE.
16	C ----------
17	C FROM ORODRAG
18	C
19	C EXPLICIT ARGUMENTS :
20	C --------------------
21	C ==== INPUTS ===
22	C ==== OUTPUTS ===
23	C
24	C IMPLICIT ARGUMENTS : NONE
25	C --------------------
26	C
27	C METHOD.
28	C -------
29	C
30	C
31	C EXTERNALS.
32	C ----------
33	C
34	C
35	C REFERENCE.
36	C ----------
37	C
38	C SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
39	C
40	C AUTHOR.
41	C -------
42	C
43	C MODIFICATIONS.
44	C --------------
45	C F.LOTT FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993
46	C
47	C-----------------------------------------------------------------------
48	use dimradmars_mod, only: ndlo2
49	USE comcstfi_h
50	implicit none
51	C
52
53	integer klon,klev,kidia,kfdia
54
55	#include "yoegwd.h"
56
57	C-----------------------------------------------------------------------
58	C
59	C* 0.1 ARGUMENTS
60	C ---------
61	C
62	INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2),KSECT(NDLO2),
63	* KKHLIM(NDLO2),KTEST(NDLO2),kkenvh(NDLO2)
64
65	C
66	REAL PAPHM1(NDLO2,KLEV+1),PAPM1(NDLO2,KLEV),PUM1(NDLO2,KLEV),
67	* PVM1(NDLO2,KLEV),PTM1(NDLO2,KLEV),PGEOM1(NDLO2,KLEV),
68	* PRHO(NDLO2,KLEV+1),PRI(NDLO2,KLEV+1),PSTAB(NDLO2,KLEV+1),
69	* PTAU(NDLO2,KLEV+1),PVPH(NDLO2,KLEV+1),ppsi(NDLO2,klev+1),
70	* pzdep(NDLO2,klev)
71	REAL PULOW(NDLO2),PVLOW(NDLO2),ptheta(NDLO2),pgamma(NDLO2),
72	* pnu(NDLO2),
73	* pd1(NDLO2),pd2(NDLO2),pdmod(NDLO2)
74	real pvaror(NDLO2)
75	C
76	C-----------------------------------------------------------------------
77	C
78	C* 0.2 LOCAL ARRAYS
79	C ------------
80	C
81	C
82	LOGICAL LL1(NDLO2,klev+1)
83	integer kknu(NDLO2),kknu2(NDLO2),kknub(NDLO2),kknul(NDLO2),
84	* kentp(NDLO2),ncount(NDLO2)
85	C
86	REAL ZHCRIT(NDLO2,klev),ZNCRIT(NDLO2,klev),
87	* ZVPF(NDLO2,klev), ZDP(NDLO2,klev)
88	REAL ZNORM(NDLO2),zpsi(NDLO2),zb(NDLO2),zc(NDLO2),
89	* zulow(NDLO2),zvlow(NDLO2),znup(NDLO2),znum(NDLO2)
90	C
91	c declarations pour "implicit none"
92	integer jk,jl,ilevm1,ilevm2,ilevh
93	real zu,zphi,zcons1,zcons2,zcons3,zwind,zdwind,zhgeo
94	real zvt1,zvt2,zst,zvar,zdelp,zstabm,zstabp,zrhom,zrhop
95	real alpha,zggeenv,zggeom1,zgvar
96	logical lo
97
98	C ------------------------------------------------------------------
99	C
100	C* 1. INITIALIZATION
101	C --------------
102	C
103	c print *,' entree gwsetup'
104	100 CONTINUE
105	C
106	C ------------------------------------------------------------------
107	C
108	C* 1.1 COMPUTATIONAL CONSTANTS
109	C -----------------------
110	C
111
112	kidia=1
113	kfdia=klon
114
115	110 CONTINUE
116	C
117	ILEVM1=KLEV-1
118	ILEVM2=KLEV-2
119	ILEVH =KLEV/3
120	C
121	ZCONS1=1./r
122	cold ZCONS2=G**2/CPD
123	ZCONS2=g**2/cpp
124	cold ZCONS3=1.5*API
125	ZCONS3=1.5*PI
126	C
127	C
128	C ------------------------------------------------------------------
129	C
130	C* 2.
131	C --------------
132	C
133	200 CONTINUE
134	C
135	C ------------------------------------------------------------------
136	C
137	C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
138	C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
139	C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
140	C
141	C
142	C
143	DO 2001 JL=kidia,kfdia
144	kknu(JL) =klev
145	kknu2(JL) =klev
146	kknub(JL) =klev
147	kknul(JL) =klev
148	pgamma(JL) =max(pgamma(jl),gtsec)
149	ll1(jl,klev+1)=.false.
150	2001 CONTINUE
151	C
152	C* DEFINE TOP OF LOW LEVEL FLOW
153	C ----------------------------
154	DO 2002 JK=KLEV,ilevh,-1
155	DO 2003 JL=kidia,kfdia
156	LO=(PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)).GE.GSIGCR
157	IF(LO) THEN
158	KKCRIT(JL)=JK
159	ENDIF
160	ZHCRIT(JL,JK)=4.*pvaror(JL)
161	ZHGEO=PGEOM1(JL,JK)/g
162	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
163	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
164	kknu(JL)=JK
165	ENDIF
166	2003 CONTINUE
167	2002 CONTINUE
168	DO 2004 JK=KLEV,ilevh,-1
169	DO 2005 JL=kidia,kfdia
170	ZHCRIT(JL,JK)=3.*pvaror(JL)
171	ZHGEO=PGEOM1(JL,JK)/g
172	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
173	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
174	kknu2(JL)=JK
175	ENDIF
176	2005 CONTINUE
177	2004 CONTINUE
178	DO 2006 JK=KLEV,ilevh,-1
179	DO 2007 JL=kidia,kfdia
180	ZHCRIT(JL,JK)=2.*pvaror(JL)
181	ZHGEO=PGEOM1(JL,JK)/g
182	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
183	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
184	kknub(JL)=JK
185	ENDIF
186	2007 CONTINUE
187	2006 CONTINUE
188	DO 2008 JK=KLEV,ilevh,-1
189	DO 2009 JL=kidia,kfdia
190	ZHCRIT(JL,JK)=pvaror(JL)
191	ZHGEO=PGEOM1(JL,JK)/g
192	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
193	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
194	kknul(JL)=JK
195	ENDIF
196	2009 CONTINUE
197	2008 CONTINUE
198	C
199	do 2010 jl=kidia,kfdia
200	kknu(jl)=min(kknu(jl),nktopg)
201	kknub(jl)=min(kknub(jl),nktopg)
202	if(kknub(jl).eq.nktopg) kknul(jl)=klev
203	C
204	C CHANGE IN HERE TO STOP KKNUL=KKNUB
205	C
206	if(kknul(jl).le.kknub(jl)) kknul(jl)=nktopg
207	2010 continue
208	C
209
210	210 CONTINUE
211	C
212	C
213	CC* INITIALIZE VARIOUS ARRAYS
214	C
215	DO 2107 JL=kidia,kfdia
216	PRHO(JL,KLEV+1) =0.0
217	PSTAB(JL,KLEV+1) =0.0
218	PSTAB(JL,1) =0.0
219	PRI(JL,KLEV+1) =9999.0
220	ppsi(JL,KLEV+1) =0.0
221	PRI(JL,1) =0.0
222	PVPH(JL,1) =0.0
223	PULOW(JL) =0.0
224	PVLOW(JL) =0.0
225	zulow(JL) =0.0
226	zvlow(JL) =0.0
227	KKCRITH(JL) =KLEV
228	KKenvH(JL) =KLEV
229	Kentp(JL) =KLEV
230	KCRIT(JL) =1
231	ncount(JL) =0
232	ll1(JL,klev+1) =.false.
233	2107 CONTINUE
234	C
235	C* DEFINE LOW-LEVEL FLOW
236	C ---------------------
237	C
238	DO 223 JK=KLEV,2,-1
239	DO 222 JL=kidia,kfdia
240	IF(KTEST(JL).EQ.1) THEN
241	ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
242	PRHO(JL,JK)=2.PAPHM1(JL,JK)ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
243	PSTAB(JL,JK)=2.ZCONS2/(PTM1(JL,JK)+PTM1(JL,JK-1))
244	* (1.-cppPRHO(JL,JK)(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
245	PSTAB(JL,JK)=MAX(PSTAB(JL,JK),GSSEC)
246	ENDIF
247	222 CONTINUE
248	223 CONTINUE
249	C
250	C********************************************************************
251	C
252	C* DEFINE blocked FLOW
253	C -------------------
254	DO 2115 JK=klev,ilevh,-1
255	DO 2116 JL=kidia,kfdia
256	if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
257	pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
258	pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
259	end if
260	2116 CONTINUE
261	2115 CONTINUE
262	DO 2110 JL=kidia,kfdia
263	pulow(JL)=pulow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
264	pvlow(JL)=pvlow(JL)/(PAPHM1(JL,Kknul(jl)+1)-PAPHM1(JL,kknub(jl)))
265	ZNORM(JL)=MAX(SQRT(PULOW(JL)2+PVLOW(JL)2),GVSEC)
266	PVPH(JL,KLEV+1)=ZNORM(JL)
267	2110 CONTINUE
268	C
269	C***** SETUP OROGRAPHY AXES AND DEFINE PLANE OF PROFILES *****
270	C
271	DO 2112 JL=kidia,kfdia
272	LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
273	IF(LO) THEN
274	ZU=PULOW(JL)+2.*GVSEC
275	ELSE
276	ZU=PULOW(JL)
277	ENDIF
278	Zphi=ATAN(PVLOW(JL)/ZU)
279	ppsi(jl,klev+1)=ptheta(jl)*pi/180.-zphi
280	zb(jl)=1.-0.18pgamma(jl)-0.04pgamma(jl)**2
281	zc(jl)=0.48pgamma(jl)+0.3pgamma(jl)**2
282	pd1(jl)=zb(jl)-(zb(jl)-zc(jl))(sin(ppsi(jl,klev+1))*2)
283	pd2(jl)=(zb(jl)-zc(jl))sin(ppsi(jl,klev+1))cos(ppsi(jl,klev+1))
284	pdmod(jl)=sqrt(pd1(jl)2+pd2(jl)2)
285	2112 CONTINUE
286	C
287	C ********** DEFINE FLOW IN PLANE OF LOWLEVEL STRESS ***********
288	C
289	DO 213 JK=1,KLEV
290	DO 212 JL=kidia,kfdia
291	IF(KTEST(JL).EQ.1) THEN
292	ZVt1 =PULOW(JL)PUM1(JL,JK)+PVLOW(JL)PVM1(JL,JK)
293	ZVt2 =-PvLOW(JL)PUM1(JL,JK)+PuLOW(JL)PVM1(JL,JK)
294	ZVPF(JL,JK)=(zvt1pd1(jl)+zvt2pd2(JL))/(znorm(jl)*pdmod(jl))
295	ENDIF
296	PTAU(JL,JK) =0.0
297	Pzdep(JL,JK) =0.0
298	Ppsi(JL,JK) =0.0
299	ll1(JL,JK) =.FALSE.
300	212 CONTINUE
301	213 CONTINUE
302	DO 215 JK=2,KLEV
303	DO 214 JL=kidia,kfdia
304	IF(KTEST(JL).EQ.1) THEN
305	ZDP(JL,JK)=PAPM1(JL,JK)-PAPM1(JL,JK-1)
306	PVPH(JL,JK)=((PAPHM1(JL,JK)-PAPM1(JL,JK-1))*ZVPF(JL,JK)+
307	* (PAPM1(JL,JK)-PAPHM1(JL,JK))*ZVPF(JL,JK-1))
308	* /ZDP(JL,JK)
309	IF(PVPH(JL,JK).LT.GVSEC) THEN
310	PVPH(JL,JK)=GVSEC
311	KCRIT(JL)=JK
312	ENDIF
313	ENDIF
314	214 CONTINUE
315	215 CONTINUE
316	C
317	C
318	C* 2.2 BRUNT-VAISALA FREQUENCY AND DENSITY AT HALF LEVELS.
319	C
320	220 CONTINUE
321	C
322	DO 2211 JK=ilevh,KLEV
323	DO 221 JL=kidia,kfdia
324	IF(KTEST(JL).EQ.1) THEN
325	IF(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) THEN
326	ZST=ZCONS2/PTM1(JL,JK)(1.-cppPRHO(JL,JK)*
327	* (PTM1(JL,JK)-PTM1(JL,JK-1))/ZDP(JL,JK))
328	PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)+ZST*ZDP(JL,JK)
329	PSTAB(JL,KLEV+1)=MAX(PSTAB(JL,KLEV+1),GSSEC)
330	PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)+PAPHM1(JL,JK)2.ZDP(JL,JK)
331	* *ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
332	ENDIF
333	ENDIF
334	221 CONTINUE
335	2211 CONTINUE
336	C
337	DO 2212 JL=kidia,kfdia
338	C*****************************************************************************
339	C
340	C O.K. THERE IS A POSSIBLE PROBLEM HERE. IF KKNUL=KKNUB THEN
341	C DIVISION BY ZERO OCCURS. I HAVE PUT A FIX IN HERE BUT WILL ASK FRANCOIS
342	C LOTT ABOUT IT IN PARIS.
343	C
344	C MAT COLLINS 30.1.96
345	C
346	C ALSO IF THIS IS THE CASE PSTAB AND PRHO ARE NOT DEFINED AT KLEV+1
347	C SO I HAVE ADDED THE ELSE
348	C
349	C*****************************************************************************
350	IF (KKNUL(JL).NE.KKNUB(JL)) THEN
351	PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
352	* -PAPM1(JL,kknub(jl)))
353	PRHO(JL,KLEV+1)=PRHO(JL,KLEV+1)/(PAPM1(JL,Kknul(jl))
354	* -PAPM1(JL,kknub(jl)))
355	ELSE
356	WRITE(,) 'OROSETUP: KKNUB=KKNUL= ',KKNUB(JL),' AT JL= ',JL
357	PSTAB(JL,KLEV+1)=PSTAB(JL,KLEV)
358	PRHO(JL,KLEV+1)=PRHO(JL,KLEV)
359	ENDIF
360	ZVAR=PVARor(JL)
361	2212 CONTINUE
362	C
363	C* 2.3 MEAN FLOW RICHARDSON NUMBER.
364	C* AND CRITICAL HEIGHT FOR FROUDE LAYER
365	C
366	230 CONTINUE
367	C
368	DO 232 JK=2,KLEV
369	DO 231 JL=kidia,kfdia
370	IF(KTEST(JL).EQ.1) THEN
371	ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
372	PRI(JL,JK)=PSTAB(JL,JK)*(ZDP(JL,JK)
373	* /(gPRHO(JL,JK)ZDWIND))**2
374	PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
375	ENDIF
376	231 CONTINUE
377	232 CONTINUE
378	C
379	C
380	C* DEFINE TOP OF 'envelope' layer
381	C ----------------------------
382
383	DO 233 JL=kidia,kfdia
384	pnu (jl)=0.0
385	znum(jl)=0.0
386	233 CONTINUE
387
388	DO 234 JK=2,KLEV-1
389	DO 234 JL=kidia,kfdia
390
391	IF(KTEST(JL).EQ.1) THEN
392
393	IF (JK.GE.KKNU2(JL)) THEN
394
395	ZNUM(JL)=PNU(JL)
396	ZWIND=(pulow(JL)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
397	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
398	ZWIND=max(sqrt(zwind**2),gvsec)
399	ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
400	ZSTABM=SQRT(MAX(PSTAB(JL,JK ),GSSEC))
401	ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
402	ZRHOM=PRHO(JL,JK )
403	ZRHOP=PRHO(JL,JK+1)
404	PNU(JL) = PNU(JL) + (ZDELP/g)*
405	* ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND
406	IF((ZNUM(JL).LE.GFRCRIT).AND.(PNU(JL).GT.GFRCRIT)
407	* .AND.(KKENVH(JL).EQ.KLEV))
408	* KKENVH(JL)=JK
409
410	ENDIF
411
412	ENDIF
413
414	234 CONTINUE
415
416	C CALCULATION OF A DYNAMICAL MIXING HEIGHT FOR THE BREAKING
417	C OF GRAVITY WAVES:
418
419
420	DO 235 JL=kidia,kfdia
421	znup(jl)=0.0
422	znum(jl)=0.0
423	235 CONTINUE
424
425	DO 236 JK=KLEV-1,2,-1
426	DO 236 JL=kidia,kfdia
427
428	IF(KTEST(JL).EQ.1) THEN
429
430	IF (JK.LT.KKENVH(JL)) THEN
431
432	ZNUM(JL)=ZNUP(JL)
433	ZWIND=(pulow(JL)pum1(jl,jk)+pvlow(jl)pvm1(jl,jk))/
434	* max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
435	ZWIND=max(sqrt(zwind**2),gvsec)
436	ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
437	ZSTABM=SQRT(MAX(PSTAB(JL,JK ),GSSEC))
438	ZSTABP=SQRT(MAX(PSTAB(JL,JK+1),GSSEC))
439	ZRHOM=PRHO(JL,JK )
440	ZRHOP=PRHO(JL,JK+1)
441	ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*
442	* ((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND
443	IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5)
444	* .AND.(KKCRITH(JL).EQ.KLEV))
445	* KKCRITH(JL)=JK
446
447	ENDIF
448
449	ENDIF
450
451	236 CONTINUE
452
453	DO 237 JL=KIDIA,KFDIA
454	KKCRITH(JL)=MIN0(KKCRITH(JL),KKNU(JL))
455	237 CONTINUE
456	c
457	c directional info for flow blocking *************************
458	c
459	do 251 jk=ilevh,klev
460	DO 252 JL=kidia,kfdia
461	IF(jk.ge.kkenvh(jl)) THEN
462	LO=(PUm1(JL,jk).LT.GVSEC).AND.(PUm1(JL,jk).GE.-GVSEC)
463	IF(LO) THEN
464	ZU=PUm1(JL,jk)+2.*GVSEC
465	ELSE
466	ZU=PUm1(JL,jk)
467	ENDIF
468	Zphi=ATAN(PVm1(JL,jk)/ZU)
469	ppsi(jl,jk)=ptheta(jl)*pi/180.-zphi
470	end if
471	252 continue
472	251 CONTINUE
473	c forms the vertical 'leakiness' **************************
474
475	alpha=3.
476
477	DO 254 JK=ilevh,klev
478	DO 253 JL=kidia,kfdia
479	IF(jk.ge.kkenvh(jl)) THEN
480	zggeenv=AMAX1(1.,
481	* (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)
482	zggeom1=AMAX1(pgeom1(jl,jk),1.)
483	zgvar=amax1(pvaror(jl)*g,1.)
484	pzdep(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))
485	end if
486	253 CONTINUE
487	254 CONTINUE
488
489	260 CONTINUE
490
491	RETURN
492	END

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

Download in other formats:

Original Format