Context Navigation

calcpar.f90 @ b398fb6

FPv9.3.1FPv9.3.1b_testingFPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair

Last change on this file since b398fb6 was 496c607, checked in by Don Morton <Don.Morton@…>, 8 years ago

Initial commit of FPv9.3.1

Currently, this is a clone of snapshot FPv9.3.0

Property mode set to 100644

File size: 12.9 KB

Line
1	!**********************************************************************
2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
5	! *
6	! This file is part of FLEXPART. *
7	! *
8	! FLEXPART is free software: you can redistribute it and/or modify *
9	! it under the terms of the GNU General Public License as published by*
10	! the Free Software Foundation, either version 3 of the License, or *
11	! (at your option) any later version. *
12	! *
13	! FLEXPART is distributed in the hope that it will be useful, *
14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16	! GNU General Public License for more details. *
17	! *
18	! You should have received a copy of the GNU General Public License *
19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
20	!**********************************************************************
21
22	subroutine calcpar_ecmwf(n,uuh,vvh,pvh)
23	! i i i o
24	!*****************************************************************************
25	! *
26	! Computation of several boundary layer parameters needed for the *
27	! dispersion calculation and calculation of dry deposition velocities. *
28	! All parameters are calculated over the entire grid. *
29	! *
30	! Author: A. Stohl *
31	! *
32	! 21 May 1995 *
33	! *
34	! ------------------------------------------------------------------ *
35	! Petra Seibert, Feb 2000: *
36	! convection scheme: *
37	! new variables in call to richardson *
38	! *
39	!*****************************************************************************
40	! Changes, Bernd C. Krueger, Feb. 2001:
41	! Variables tth and qvh (on eta coordinates) in common block
42	!*****************************************************************************
43	! Changes Arnold, D. and Morton, D. (2015): *
44	! -- description of local and common variables *
45	!*****************************************************************************
46	! Functions: *
47	! scalev computation of ustar *
48	! obukhov computatio of Obukhov length *
49	! *
50	!*****************************************************************************
51
52	use par_mod
53	use com_mod
54
55	implicit none
56
57	!***********************************************************************
58	! Subroutine Parameters: *
59	! input *
60	! n temporal index for meteorological fields (1 to 3)*
61	! uuh,vvh, wwh wind components in ECMWF model levels *
62	integer :: n
63	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
64	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
65	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
66
67	!***********************************************************************
68	! Local variables *
69	! *
70	! ttlev
71	! qvlev
72	! * obukhov_ecmwf subroutine/function to calculate Obukhov length *
73	! * scalev subroutine/function to calculate ustar *
74	! ol obukhov length
75	! hmixplus maximum lifting from availiable kinetic enrgy *
76	! ulev, vlev wind speed at model levels *
77	! ew subroutine/function to calculate saturation *
78	! water vaport for a given air temperature *
79	! rh relative humidity at surface *
80	! vd deposition velocity from all species *
81	! subsceff excess hmin due to subgrid effects *
82	! ylat temporary latitude *
83	! altmin minimum height of the tropopause *
84	! tvoldm pold, zold temporary variables to keep previous values *
85	! pint pressure on model levels *
86	! tv virtual temperature on model levels *
87	! zlev height of model levels *
88	real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov_ecmwf,scalev,ol,hmixplus
89	real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
90	real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
91
92	! Other variables:
93	! ix,jy,kz,i,lz,kzmin loop control indices in each direction *
94	integer :: ix,jy,i,kz,lz,kzmin
95
96	!***********************************************************************
97
98	!***********************************************************************
99	! Local constants *
100	real,parameter :: const=r_air/ga
101	!***********************************************************************
102
103
104	!***********************************************************************
105	! Global variables *
106	! from par_mod and com_mod *
107	! ustar [m/s] friction velocity *
108	! oli [m] inverse Obukhov length (1/L) *
109	! hmix [m] mixing height *
110	! wstar [m/s] convective velocity scale *
111	! ustar [m/s] friction velocity *
112	! z0 roughness length for the landuse classes *
113	! tropopause [m] altitude of thermal tropopause *
114	! nx, ny actual dimensions of wind fields in x and y direction *
115	! dx, dy grid distances in x,y direction *
116	! akm, bkm coefficients which regulate vertical discretization of ecmwf*
117	! akz, bkz model discretization coeffizients at the centre of layers *
118	! ps surface pressure *
119	! tt2 2-m temperature *
120	! tt2d 2-m dew temperature *
121	! sshf surface sensible heat flux *
122	! surfstr surface stress *
123	! convprec convective precip *
124	! lsprec large scale precip *
125	! sd snow depth *
126	! ssr surface solar radiation *
127	! xlon0, ylat0 geographical longitude/latitude of lower left grid point*
128	!
129	!***********************************************************************
130
131	!-----------------------------------------------------------------------------
132
133
134	!write(,) 'in calcpar writting snowheight'
135	!***********************************
136	! for test: write out snow depths
137
138	! open(4,file='slandusetest',form='formatted')
139	! do 5 ix=0,nxmin1
140	!5 write (4,*) (sd(ix,jy,1,n),jy=0,nymin1)
141	! close(4)
142
143
144	! Loop over entire grid
145	!**********************
146
147	do jy=0,nymin1
148
149	! Set minimum height for tropopause
150	!**********************************
151
152	ylat=ylat0+real(jy)*dy
153	if ((ylat.ge.-20.).and.(ylat.le.20.)) then
154	altmin = 5000.
155	else
156	if ((ylat.gt.20.).and.(ylat.lt.40.)) then
157	altmin=2500.+(40.-ylat)*125.
158	else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
159	altmin=2500.+(40.+ylat)*125.
160	else
161	altmin=2500.
162	endif
163	endif
164
165	do ix=0,nxmin1
166
167	! 1) Calculation of friction velocity
168	!************************************
169
170	ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
171	td2(ix,jy,1,n),surfstr(ix,jy,1,n))
172	if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
173
174	! 2) Calculation of inverse Obukhov length scale
175	!***********************************************
176
177	ol=obukhov_ecmwf(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
178	tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm)
179	if (ol.ne.0.) then
180	oli(ix,jy,1,n)=1./ol
181	else
182	oli(ix,jy,1,n)=99999.
183	endif
184
185
186	! 3) Calculation of convective velocity scale and mixing height
187	!**************************************************************
188
189	do i=1,nuvz
190	ulev(i)=uuh(ix,jy,i)
191	vlev(i)=vvh(ix,jy,i)
192	ttlev(i)=tth(ix,jy,i,n)
193	qvlev(i)=qvh(ix,jy,i,n)
194	end do
195
196	call richardson_ecmwf(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
197	ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
198	td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus)
199
200	if(lsubgrid.eq.1) then
201	subsceff=min(excessoro(ix,jy),hmixplus)
202	else
203	subsceff=0.0
204	endif
205	!
206	! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
207	!
208	hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
209	hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
210	hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
211
212	! 4) Calculation of dry deposition velocities
213	!********************************************
214
215	if (DRYDEP) then
216	! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
217	! windspeed
218	z0(7)=0.016ustar(ix,jy,1,n)ustar(ix,jy,1,n)/ga
219
220	! Calculate relative humidity at surface
221	!***************************************
222	rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
223
224	call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
225	tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
226	ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
227	sd(ix,jy,1,n),vd)
228
229	do i=1,nspec
230	vdep(ix,jy,i,n)=vd(i)
231	end do
232
233	endif
234
235	!******************************************************
236	! Calculate height of thermal tropopause (Hoinka, 1997)
237	!******************************************************
238
239	! 1) Calculate altitudes of ECMWF model levels
240	!*********************************************
241
242	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
243	ps(ix,jy,1,n))
244	pold=ps(ix,jy,1,n)
245	zold=0.
246	do kz=2,nuvz
247	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers
248	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
249
250	if (abs(tv-tvold).gt.0.2) then
251	zlev(kz)=zold+constlog(pold/pint)(tv-tvold)/log(tv/tvold)
252	else
253	zlev(kz)=zold+constlog(pold/pint)tv
254	endif
255	tvold=tv
256	pold=pint
257	zold=zlev(kz)
258	end do
259
260	! 2) Define a minimum level kzmin, from which upward the tropopause is
261	! searched for. This is to avoid inversions in the lower troposphere
262	! to be identified as the tropopause
263	!************************************************************************
264
265	do kz=1,nuvz
266	if (zlev(kz).ge.altmin) then
267	kzmin=kz
268	goto 45
269	endif
270	end do
271	45 continue
272
273	! 3) Search for first stable layer above minimum height that fulfills the
274	! thermal tropopause criterion
275	!************************************************************************
276
277	do kz=kzmin,nuvz
278	do lz=kz+1,nuvz
279	if ((zlev(lz)-zlev(kz)).gt.2000.) then
280	if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
281	(zlev(lz)-zlev(kz))).lt.0.002) then
282	tropopause(ix,jy,1,n)=zlev(kz)
283	goto 51
284	endif
285	goto 50
286	endif
287	end do
288	50 continue
289	end do
290	51 continue
291
292
293	end do
294	end do
295
296	! Calculation of potential vorticity on 3-d grid
297	!***********************************************
298
299	call calcpv(n,uuh,vvh,pvh)
300
301
302	end subroutine calcpar_ecmwf

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

Download in other formats:

Original Format