Context Navigation

orosetup.F90 @ 2740

Last change on this file since 2740 was 2651, checked in by emillour, 3 years ago

Mars GCM:

turn sugwd.F to sugwd.F90 with extra comments
turn yoegwd.h into a module

JL+EM

File size: 20.6 KB

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

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

Download in other formats:

Original Format