Context Navigation

orosetup.F90

Last change on this file was 3754, checked in by afalco, 3 months ago
Pluto: imported orographic gravity waves from Mars. AF
File size: 19.9 KB

Line
1	MODULE orosetup_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE OROSETUP( ngrid, nlayer, ktest, pplev, pplay, pu, pv, pt, zgeom, &
8	pvar,pthe, pgam, &
9	!output in capital
10	IKCRIT, IKCRITH, ICRIT, IKENVH,IKNU,IKNU2, &
11	!ISECT, IKHLIM, not used
12	ZRHO,PRI,BV,ZTAU,ZVPH,ZPSI,ZZDEP,ZNU,ZD1,ZD2,ZDMOD, &
13	PULOW, PVLOW)
14	!---------------------------------------------------------------------------------------------
15	!!**** GWSETUP! Computes low level stresses using subcritical and super critical forms.
16	! As well as, computes anisotropy coefficient as measure of orographic two-dimensionality
17	! F.LOTT FOR THE NEW-GWDRAG SCHEME NOVEMBER 1993!
18	!--
19	! REFERENCE.
20	! 1. SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
21	! 2. Lott, F., & Miller, M. J. (1997). A new subgrid‐scale orographic drag parametrization:
22	! Its formulation and testing.Quarterly Journal of the Royal cMeteorological Society,
23	! 123(537), 101-127.
24	!--
25	! MODIFICATIONS.
26	! 1.Rewiten by J.liu 03/03/2022
27	!-----------------------------------------------------------------------
28	use dimphy, only: ndomainsz
29	use comcstfi_mod, only: cpp, g, r, pi
30	use yoegwd_h, only: gfrcrit, grcrit, gsigcr, gssec, gtsec, gvsec
31	use yoegwd_h, only: nktopg
32
33	implicit none
34
35	! 0. DECLARATIONS:
36
37	! 0.1 inputs:
38	integer,intent(in):: ngrid ! number of atmospheric columns
39	integer,intent(in):: nlayer ! number of atmospheric layers
40	INTEGER,intent(in):: ktest(ndomainsz) ! map of calling points
41
42	real, intent(in) :: pplev(ndomainsz,nlayer+1)! Pressure at 1/2 levels(Pa) (has been inversed by DRAG_NORO=inv_pplev)
43	real, intent(in) :: pplay(ndomainsz,nlayer) ! Pressure at full levels(Pa) (has been inversed by DRAG_NORO=inv_pplay)
44	real, intent(in) :: pu(ndomainsz,nlayer) ! Zonal wind at full levels(m/s) (has been inversed by DRAG_NORO, =inv_pu)
45	real, intent(in) :: pv(ndomainsz,nlayer) ! Meridional winds at full levels(m/s)(has been inversed by DRAG_NORO, =inv_pv)
46	real, intent(in) :: pt(ndomainsz,nlayer) ! Temperature at full levels(m/s) (has been inversed by DRAG_NORO=inv_pt)
47	real, intent(in) :: zgeom(ndomainsz,nlayer) ! Geopotetial height
48	real, intent(in) :: pvar(ndomainsz) ! Sub-grid scale standard deviation
49	real, intent(in) :: pthe(ndomainsz) ! Sub-grid scale principal axes angle
50	real, intent(inout) :: pgam(ndomainsz) ! Sub-grid scale anisotropy
51
52	! 0.2 outputs:
53	INTEGER,intent(out):: IKCRIT(ndomainsz) ! top of low level flow height
54	INTEGER,intent(out):: IKCRITH(ndomainsz) ! dynamical mixing height for the breaking of gravity waves
55	INTEGER,intent(out):: ICRIT(ndomainsz) ! Critical layer where orographic GW breaks
56	! INTEGER,intent(out):: ISECT(ndomainsz) ! not used
57	! INTEGER,intent(out):: IKHLIM(ndomainsz) ! not used
58	INTEGER,intent(out):: IKENVH(ndomainsz) ! Top of the blocked layer
59	INTEGER,intent(out):: IKNU(ndomainsz) ! 4*pvar layer
60	INTEGER,intent(out):: IKNU2(ndomainsz) ! 3*pvar layer
61
62	REAL, intent(out):: ZRHO(ndomainsz,nlayer+1) ! Density at 1/2 level
63	REAL, intent(out):: PRI(ndomainsz,nlayer+1) ! Mean flow richardson number at 1/2 level
64	REAL, intent(out):: BV(ndomainsz,nlayer+1) ! Brunt–Väisälä frequency at 1/2 level
65	REAL, intent(out):: ZTAU(ndomainsz,nlayer+1) ! Gravity wave stress. Set to 0.0 here and will calculate in GWSTRESS later
66	REAL, intent(out):: ZVPH(ndomainsz,nlayer+1) ! Low level wind speed U_H
67	REAL, intent(out):: ZPSI(ndomainsz,nlayer+1) ! The angle between the incident flow direction and the normal ridge direction pthe
68	REAL, intent(out):: ZZDEP(ndomainsz,nlayer) ! dp by full level
69
70	REAL, intent(out):: PULOW(ndomainsz) ! Low level zonal wind
71	REAL, intent(out):: PVLOW(ndomainsz) ! Low level meridional wind
72	REAL, intent(out):: ZNU(ndomainsz) ! A critical value see equation 9
73	REAL, intent(out):: ZD1(ndomainsz) ! Bcos^2(psi)-Csin^2(psi) see equation 17 or 18
74	REAL, intent(out):: ZD2(ndomainsz) ! (B-C)sin(psi)cos(psi) see equation 17 or 18
75	REAL, intent(out):: ZDMOD(ndomainsz) ! sqrt(zd1^2+zd2^2)
76
77	!0.3 Local arrays
78	integer IKNUb(ndomainsz) ! 2*pvar layer
79	integer IKNUl(ndomainsz) ! 1*pvar layer
80	integer kentp(ndomainsz) ! initialized value but never used
81	integer ncount(ndomainsz) ! initialized value but never used
82
83	REAL ZHCRIT(ndomainsz,nlayer) ! tag for 1pvar, 2pvar,3pvar and 4pvar, pvar is mu means SD
84	! REAL ZNCRIT(ndomainsz,nlayer) ! not used
85	REAL ZVPF(ndomainsz,nlayer) ! Flow in plane of low level stress. Seems a unit vector
86	REAL ZDP(ndomainsz,nlayer) ! dp differitial of pressure
87
88	REAL ZNORM(ndomainsz) ! The norm ridge of a moutain?
89	REAL zb(ndomainsz) ! Parameter B in eqution 17
90	REAL zc(ndomainsz) ! Parameter C in eqution 17
91	REAL zulow(ndomainsz) ! initialized value but never used
92	REAL zvlow(ndomainsz) ! initialized value but never used
93	REAL znup(ndomainsz) ! znu in top of 1/2 level
94	REAL znum(ndomainsz) ! znu in bottom of 1/2 level
95
96	integer jk,jl
97	integer ilevm1 !=nlayer-1
98	integer ilevm2 !=nlayer-2
99	integer ilevh !=nalyer/3
100	INTEGER kidia !=1
101	INTEGER kfdia !=ngrid
102	real zcons1 !=1/r
103	real zcons2 !=g^2/cpp
104	real zcons3 !=1.5*pi
105	real zphi ! direction of the incident flow
106	real zhgeo ! Height calculated by geopotential/g
107	real zu ! Low level zonal wind (to denfine the dirction of background wind)
108	real zwind,zdwind
109	real zvt1,zvt2,zst,zvar
110	real zdelp !dp differitial of pressure
111	! variables for bv and density rho
112	real zstabm,zstabp,zrhom,zrhop
113	! real alpha !=3. but never used
114	real zggeenv,zggeom1,zgvar
115	logical lo
116	LOGICAL LL1(ndomainsz,nlayer+1)
117
118	!--------------------------------------------------------------------------------
119	! 1. INITIALIZATION
120	!--------------------------------------------------------------------------------
121
122	!* 1.1 COMPUTATIONAL CONSTANTS
123	kidia=1
124	kfdia=ngrid
125	ILEVM1=nlayer-1
126	ILEVM2=nlayer-2
127	ILEVH =nlayer/3 !!!! maybe not enough for Mars, need check later
128	ZCONS1=1./r
129	ZCONS2=g**2/cpp
130	ZCONS3=1.5*PI
131
132	!------------------------------------------------------------------------------------------------------
133	! 2. Compute all the critical levels and the coeffecients of anisotropy
134	!-----------------------------------------------------------------------------------------------------
135
136	! 2.1 Define low level wind, project winds in plane of low level wind,
137	! determine sector in which to take the variance and set indicator for critical levels.
138	DO JL=kidia,kfdia
139	! initialize all the height into surface (notice the layers have been inversed by preious rountines)
140	IKNU(JL) =nlayer
141	IKNU2(JL) =nlayer
142	IKNUb(JL) =nlayer
143	IKNUl(JL) =nlayer
144	pgam(JL) =max(pgam(jl),gtsec) ! gtsec is from yoegwd.h which is a common variable
145	ll1(jl,nlayer+1)=.false.
146	end DO
147
148	! Define top of low level flow (since pressure, zonal and meridional wind have been inversed
149	! the process is to find the layer from surface (nlayer) to some levels ) by searching several
150	! altitude scope
151
152	! using 4 times sub-grid scale deviation as the threahold of the critical height
153	DO JK=nlayer,ilevh,-1 ! ilevh=nlayer/3=16
154	DO JL=kidia,kfdia ! jl=1:ngrid
155	! To found the layer of the "top of low level flow"
156	LO=(pplev(JL,JK)/pplev(JL,nlayer+1)).GE.GSIGCR
157	IF(LO) THEN
158	IKCRIT(JL)=JK
159	ENDIF
160	ZHCRIT(JL,JK)=4.*pvar(JL) !
161	! use geopotetial denfination to get geoheight[in meters] of the layer
162	ZHGEO=zgeom(JL,JK)/g
163	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
164	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
165	IKNU(JL)=JK
166	ENDIF
167	ENDDO
168	end DO
169
170	! using 3 times sub-grid scale deviation as the threahold of the critical height
171	DO JK=nlayer,ilevh,-1
172	DO JL=kidia,kfdia
173	ZHCRIT(JL,JK)=3.*pvar(JL)
174	ZHGEO=zgeom(JL,JK)/g
175	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
176	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
177	IKNU2(JL)=JK
178	ENDIF
179	ENDDO
180	end DO
181
182	! using 2 times sub-grid scale deviation as the threahold of the critical height
183	DO JK=nlayer,ilevh,-1
184	DO JL=kidia,kfdia
185	ZHCRIT(JL,JK)=2.*pvar(JL)
186	ZHGEO=zgeom(JL,JK)/g
187	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
188	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
189	IKNUb(JL)=JK
190	ENDIF
191	ENDDO
192	end DO
193
194	! using 1 times sub-grid scale deviation as the threahold of the critical height
195	DO JK=nlayer,ilevh,-1
196	DO JL=kidia,kfdia
197	ZHCRIT(JL,JK)=pvar(JL)
198	ZHGEO=zgeom(JL,JK)/g
199	ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
200	IF(ll1(JL,JK).NEQV.ll1(JL,JK+1)) THEN
201	IKNUl(JL)=JK
202	ENDIF
203	ENDDO
204	end DO
205	! loop to relocate the critical height to make sure everything is okay if theses
206	! levels hit the model surface or top.
207	do jl=kidia,kfdia
208	IKNU(jl)=min(IKNU(jl),nktopg) ! nktopg is a common variable from yoegwd.h
209	IKNUb(jl)=min(IKNUb(jl),nktopg)
210	if(IKNUb(jl).eq.nktopg) IKNUl(jl)=nlayer
211	! Change in here to stop IKNUl=IKNUb
212	if(IKNUl(jl).le.IKNUb(jl)) IKNUl(jl)=nktopg
213	enddo
214
215	! Initialize various arrays for the following computes
216	DO JL=kidia,kfdia
217	ZRHO(JL,nlayer+1) =0.0
218	BV(JL,nlayer+1) =0.0
219	BV(JL,1) =0.0
220	PRI(JL,nlayer+1) =9999.0
221	ZPSI(JL,nlayer+1) =0.0
222	PRI(JL,1) =0.0
223	ZVPH(JL,1) =0.0
224	PULOW(JL) =0.0
225	PVLOW(JL) =0.0
226	zulow(JL) =0.0
227	zvlow(JL) =0.0
228	IKCRITH(JL) =nlayer ! surface
229	IKENVH(JL) =nlayer ! surface
230	Kentp(JL) =nlayer ! surface
231	ICRIT(JL) =1 ! topmost layer
232	ncount(JL) =0
233	ll1(JL,nlayer+1) =.false.
234	ENDDO
235
236	! Define low-level flow Brunt–Väisälä frequency N^2, density ZRHO
237	! The incident flow passes over the mean orography is evaluated by averaging the wind,
238	! the Brunt–Väisälä frequency, and the fluid density between 1pvar and 2pvar over the
239	! model mean orography
240	DO JK=nlayer,2,-1 ! from surface to topmost-1 layer
241	DO JL=kidia,kfdia
242	IF(ktest(JL).EQ.1) THEN ! if the map of the calling points is true
243	! calcalate density and BV
244	ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1) !dp>0
245	ZRHO(JL,JK)=2.pplev(JL,JK)ZCONS1/(pt(JL,JK)+pt(JL,JK-1)) !rho=p/(r*T)
246	! Brunt–Väisälä frequency N^2. This equation for BV is illness since
247	! too many variables are used. Use N^2=g/T[1/(cpp*T)+dT/dz] to replace in the future
248	BV(JL,JK)=2.ZCONS2/(pt(JL,JK)+pt(JL,JK-1))(1.-cppZRHO(JL,JK)(pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK))
249	BV(JL,JK)=MAX(BV(JL,JK),GSSEC)
250	ENDIF
251	ENDDO
252	end DO
253
254	DO JK=nlayer,ilevh,-1
255	DO JL=kidia,kfdia
256	if(jk.ge.IKNUb(jl).and.jk.le.IKNUl(jl)) then ! IF the layer between 1pvar and 2pvar
257	! calculate the low level wind U_H
258	! pulow/pvlow at a speicfic location equals to sum of u*dp of all levels
259	! notice here dp is already a positive number
260	pulow(JL)=pulow(JL)+pu(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK))
261	pvlow(JL)=pvlow(JL)+pv(JL,JK)*(pplev(JL,JK+1)-pplev(JL,JK))
262	end if
263	ENDDO
264	end DO
265	! averaging the wind
266	DO JL=kidia,kfdia
267	! by divide dp [p differ between iknul and uknub level]
268	pulow(JL)=pulow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl)))
269	pvlow(JL)=pvlow(JL)/(pplev(JL,IKNUl(jl)+1)-pplev(JL,IKNUb(jl)))
270	! average U to get background U?
271	ZNORM(JL)=MAX(SQRT(PULOW(JL)2+PVLOW(JL)2),GVSEC)
272	ZVPH(JL,nlayer+1)=ZNORM(JL) ! The wind below the surface level (e.g., start of the 1/2 level)
273	ENDDO
274
275	! The gravity wave drag caused by the flow passes over an single elliptic mountain can be calculated
276	! by equation 17 and 18
277	DO JL=kidia,kfdia
278	LO=(PULOW(JL).LT.GVSEC).AND.(PULOW(JL).GE.-GVSEC)
279	IF(LO) THEN
280	ZU=PULOW(JL)+2.*GVSEC
281	ELSE
282	ZU=PULOW(JL)
283	ENDIF
284	! Here all physics for equation 17 and 18
285	! Direction of the incident flow
286	Zphi=ATAN(PVLOW(JL)/ZU)
287	! The angle between the incident flow direction and the normal ridge direction pthe
288	ZPSI(jl,nlayer+1)=pthe(jl)*pi/180.-zphi
289	! equation(17) parameter B and C
290	zb(jl)=1.-0.18pgam(jl)-0.04pgam(jl)**2
291	zc(jl)=0.48pgam(jl)+0.3pgam(jl)**2
292	! Bcos^2(psi)-Csin^2(psi)
293	ZD1(jl)=zb(jl)-(zb(jl)-zc(jl))(sin(ZPSI(jl,nlayer+1))*2)
294	! (B-C)sin(psi)cos(psi)
295	ZD2(jl)=(zb(jl)-zc(jl))sin(ZPSI(jl,nlayer+1))cos(ZPSI(jl,nlayer+1))
296	! squre root of tao1 and tao2 without the constant see equation 17 or 18
297	ZDMOD(jl)=sqrt(ZD1(jl)2+ZD2(jl)2)
298	ENDDO
299
300	! Define blocked flow
301	! Setup orogrphy axes and define plane of profiles
302	! Define blocked flow in plane of the low level stress
303	DO JK=1,nlayer
304	DO JL=kidia,kfdia
305	IF(ktest(JL).EQ.1) THEN
306	ZVt1 =PULOW(JL)pu(JL,JK)+PVLOW(JL)pv(JL,JK)
307	ZVt2 =-PvLOW(JL)pu(JL,JK)+PuLOW(JL)pv(JL,JK)
308	! zvpf is a normalized variable
309	ZVPF(JL,JK)=(zvt1ZD1(jl)+zvt2ZD2(JL))/(znorm(jl)*ZDMOD(jl))
310	ENDIF
311	ZTAU(JL,JK) =0.0
312	ZZDEP(JL,JK) =0.0
313	ZPSI(JL,JK) =0.0
314	ll1(JL,JK) =.FALSE.
315	ENDDO
316	end DO
317
318	DO JK=2,nlayer
319	DO JL=kidia,kfdia
320	IF(ktest(JL).EQ.1) THEN
321	ZDP(JL,JK)=pplay(JL,JK)-pplay(JL,JK-1) ! dp
322	! zvph is the U_H in equation 17 e.g. low level wind speed
323	ZVPH(JL,JK)=((pplev(JL,JK)-pplay(JL,JK-1))*ZVPF(JL,JK)+ &
324	(pplay(JL,JK)-pplev(JL,JK))*ZVPF(JL,JK-1))/ZDP(JL,JK)
325	IF(ZVPH(JL,JK).LT.GVSEC) THEN
326	ZVPH(JL,JK)=GVSEC
327	ICRIT(JL)=JK
328	ENDIF
329	endIF
330	ENDDO
331	end DO
332
333	! 2.2 Brunt-vaisala frequency and density at half levels
334
335	DO JK=ilevh,nlayer
336	DO JL=kidia,kfdia
337	IF(ktest(JL).EQ.1) THEN
338	IF(jk.ge.(IKNUb(jl)+1).and.jk.le.IKNUl(jl)) THEN
339	ZST=ZCONS2/pt(JL,JK)(1.-cppZRHO(JL,JK)* &
340	(pt(JL,JK)-pt(JL,JK-1))/ZDP(JL,JK))
341	BV(JL,nlayer+1)=BV(JL,nlayer+1)+ZST*ZDP(JL,JK)
342	BV(JL,nlayer+1)=MAX(BV(JL,nlayer+1),GSSEC)
343	ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)+pplev(JL,JK)2.ZDP(JL,JK) &
344	*ZCONS1/(pt(JL,JK)+pt(JL,JK-1))
345	ENDIF
346	endIF
347	ENDDO
348	end DO
349
350	DO JL=kidia,kfdia
351	!*****************************************************************************
352	! Okay. There is a possible problem here. If IKNUl=IKNUb then division by zero
353	! occurs. I have put a fix in here but will ask Francois lott about it in Paris.
354	! Also if this is the case BV and ZRHO are not defined at nlayer+1 so I have
355	! added the else.
356	! by: MAT COLLINS 30.1.96
357	!*****************************************************************************
358	IF (IKNUL(JL).NE.IKNUB(JL)) THEN
359	BV(JL,nlayer+1)=BV(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl)))
360	ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer+1)/(pplay(JL,IKNUl(jl))-pplay(JL,IKNUb(jl)))
361	ELSE
362	WRITE(,) 'OROSETUP: IKNUB=IKNUL= ',IKNUB(JL),' AT JL= ',JL
363	BV(JL,nlayer+1)=BV(JL,nlayer)
364	ZRHO(JL,nlayer+1)=ZRHO(JL,nlayer)
365	ENDIF
366	ZVAR=pvar(JL)
367	ENDDO !JL=kidia,kfdia
368
369	! 2.3 Mean flow richardson number and critical height for proude layer
370
371	DO JK=2,nlayer
372	DO JL=kidia,kfdia
373	IF(ktest(JL).EQ.1) THEN
374	! du
375	ZDWIND=MAX(ABS(ZVPF(JL,JK)-ZVPF(JL,JK-1)),GVSEC)
376	! Mean flow Richardson number Ri=g/rho[drho/dz / (du/dz)^2]
377	! Here dp maybe dp^2 ? Need ask Francios lott later
378	PRI(JL,JK)=BV(JL,JK)(ZDP(JL,JK)/(gZRHO(JL,JK)ZDWIND))*2
379	PRI(JL,JK)=MAX(PRI(JL,JK),GRCRIT)
380	ENDIF
381	ENDDO
382	end DO
383
384	!* Define top of 'envelope' layer
385	DO JL=kidia,kfdia
386	ZNU (jl)=0.0
387	znum(jl)=0.0
388	ENDDO
389
390	DO JK=2,nlayer-1
391	DO JL=kidia,kfdia
392	IF(ktest(JL).EQ.1) THEN
393	IF (JK.GE.IKNU2(JL)) THEN ! level lower than 3*par
394	! all codes here is to calculate equation 9
395	ZNUM(JL)=ZNU(JL)
396	ZWIND=(pulow(JL)pu(jl,jk)+pvlow(jl)pv(jl,jk))/max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
397	ZWIND=max(sqrt(zwind**2),gvsec)
398	ZDELP=pplev(JL,JK+1)-pplev(JL,JK) ! dp
399	ZSTABM=SQRT(MAX(BV(JL,JK ),GSSEC))
400	ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC))
401	ZRHOM=ZRHO(JL,JK )
402	ZRHOP=ZRHO(JL,JK+1)
403	! Equation 9. znu is a critical value to find the blocking layer
404	ZNU(JL) = ZNU(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND
405	! Found the moutain top
406	IF((ZNUM(JL).LE.GFRCRIT).AND.(ZNU(JL).GT.GFRCRIT).AND.(IKENVH(JL).EQ.nlayer)) THEN
407	IKENVH(JL)=JK
408	ENDIF
409	ENDIF ! (JK.GE.IKNU2(JL))
410	ENDIF !(ktest(JL).EQ.1)
411	ENDDO
412	endDO
413
414	! Calculation of a dynamical mixing height for the breaking of gravity waves
415	DO JL=kidia,kfdia
416	znup(jl)=0.0
417	znum(jl)=0.0
418	ENDDO
419
420	DO JK=nlayer-1,2,-1
421	DO JL=kidia,kfdia
422	IF(ktest(JL).EQ.1) THEN
423	IF (JK.LT.IKENVH(JL)) THEN
424	ZNUM(JL)=ZNUP(JL)
425	ZWIND=(pulow(JL)pu(jl,jk)+pvlow(jl)pv(jl,jk))/max(sqrt(pulow(jl)2+pvlow(jl)2),gvsec)
426	ZWIND=max(sqrt(zwind**2),gvsec)
427	ZDELP=pplev(JL,JK+1)-pplev(JL,JK)
428	ZSTABM=SQRT(MAX(BV(JL,JK ),GSSEC))
429	ZSTABP=SQRT(MAX(BV(JL,JK+1),GSSEC))
430	ZRHOM=ZRHO(JL,JK )
431	ZRHOP=ZRHO(JL,JK+1)
432	ZNUP(JL) = ZNUP(JL) + (ZDELP/g)*((zstabp/zrhop+zstabm/zrhom)/2.)/ZWIND
433	! dynamical mixing height for the breaking of gravity waves
434	IF((ZNUM(JL).LE.1.5).AND.(ZNUP(JL).GT.1.5).AND.(IKCRITH(JL).EQ.nlayer)) THEN
435	IKCRITH(JL)=JK
436	ENDIF
437	ENDIF ! (JK.LT.IKENVH(JL))
438	ENDIF ! (ktest(JL).EQ.1)
439	ENDDO
440	end DO
441
442	DO JL=KIDIA,KFDIA
443	IKCRITH(JL)=MIN0(IKCRITH(JL),IKNU(JL))
444	ENDDO
445
446	! directional info for flow blocking *************************
447	DO jk=ilevh,nlayer
448	DO JL=kidia,kfdia
449	IF(jk.ge.IKENVH(jl)) THEN
450	LO=(pu(JL,jk).LT.GVSEC).AND.(pu(JL,jk).GE.-GVSEC)
451	IF(LO) THEN
452	ZU=pu(JL,jk)+2.*GVSEC
453	ELSE
454	ZU=pu(JL,jk)
455	ENDIF
456	Zphi=ATAN(pv(JL,jk)/ZU)
457	ZPSI(jl,jk)=pthe(jl)*pi/180.-zphi
458	end IF
459	ENDDO
460	end DO
461
462	! forms the vertical 'leakiness' **************************
463	DO JK=ilevh,nlayer
464	DO JL=kidia,kfdia
465	IF(jk.ge.IKENVH(jl)) THEN
466	zggeenv=AMAX1(1.,(zgeom(jl,IKENVH(jl))+zgeom(jl,IKENVH(jl)-1))/2.)
467	zggeom1=AMAX1(zgeom(jl,jk),1.)
468	zgvar=amax1(pvar(jl)*g,1.)
469	ZZDEP(jl,jk)=SQRT((zggeenv-zggeom1)/(zggeom1+zgvar))
470	endIF
471	ENDDO
472	end DO
473
474	END SUBROUTINE OROSETUP
475
476	END MODULE orosetup_mod

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

Download in other formats:

Original Format