Context Navigation

calcpar_gfs.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: 13.3 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_gfs(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	! *
44	! CHANGE 17/11/2005 Caroline Forster NCEP GFS version *
45	! *
46	!*****************************************************************************
47	! Changes Arnold, D. and Morton, D. (2015): *
48	! -- description of local and common variables *
49	!*****************************************************************************
50	! Functions: *
51	! scalev computation of ustar *
52	! obukhov computatio of Obukhov length *
53	! *
54	!*****************************************************************************
55
56	use par_mod
57	use com_mod
58
59	implicit none
60
61	!***********************************************************************
62	! Subroutine Parameters: *
63	! input *
64	! n temporal index for meteorological fields (1 to 3)*
65	! uuh,vvh, wwh wind components in ECMWF model levels *
66	integer :: n
67	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
68	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
69	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
70
71	!***********************************************************************
72	! Local variables *
73	! *
74	! ttlev
75	! qvlev
76	! * obukhov_gfs subroutine/function to calculate Obukhov length *
77	! * scalev subroutine/function to calculate ustar *
78	! ol obukhov length
79	! hmixplus maximum lifting from availiable kinetic enrgy *
80	! ulev, vlev wind speed at model levels *
81	! ew subroutine/function to calculate saturation *
82	! water vaport for a given air temperature *
83	! rh relative humidity at surface *
84	! vd deposition velocity from all species *
85	! subsceff excess hmin due to subgrid effects *
86	! ylat temporary latitude *
87	! altmin minimum height of the tropopause *
88	! tvoldm pold, zold temporary variables to keep previous values *
89	! pint pressure on model levels *
90	! tv virtual temperature on model levels *
91	! zlev height of model levels *
92
93	real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov_gfs,scalev,ol,hmixplus
94	real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
95	real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy
96
97
98	! Other variables:
99	! ix,jy,kz,i,lz,kzmin loop control indices in each direction *
100	integer :: ix,jy,i,kz,lz,kzmin,llev
101
102	!***********************************************************************
103
104	!***********************************************************************
105	! Local constants *
106	real,parameter :: const=r_air/ga
107	!***********************************************************************
108
109
110	!***********************************************************************
111	! Global variables *
112	! from par_mod and com_mod *
113	! ustar [m/s] friction velocity *
114	! oli [m] inverse Obukhov length (1/L) *
115	! hmix [m] mixing height *
116	! wstar [m/s] convective velocity scale *
117	! ustar [m/s] friction velocity *
118	! z0 roughness length for the landuse classes *
119	! tropopause [m] altitude of thermal tropopause *
120	! nx, ny actual dimensions of wind fields in x and y direction *
121	! dx, dy grid distances in x,y direction *
122	! akm, bkm coefficients which regulate vertical discretization of ecmwf*
123	! akz, bkz model discretization coeffizients at the centre of layers *
124	! ps surface pressure *
125	! tt2 2-m temperature *
126	! tt2d 2-m dew temperature *
127	! sshf surface sensible heat flux *
128	! surfstr surface stress *
129	! convprec convective precip *
130	! lsprec large scale precip *
131	! sd snow depth *
132	! ssr surface solar radiation *
133	! xlon0, ylat0 geographical longitude/latitude of lower left grid point*
134	!
135	!***********************************************************************
136
137	!-----------------------------------------------------------------------------
138
139
140	! Loop over entire grid
141	!**********************
142
143	do jy=0,nymin1
144
145	! Set minimum height for tropopause
146	!**********************************
147
148	ylat=ylat0+real(jy)*dy
149	if ((ylat.ge.-20.).and.(ylat.le.20.)) then
150	altmin = 5000.
151	else
152	if ((ylat.gt.20.).and.(ylat.lt.40.)) then
153	altmin=2500.+(40.-ylat)*125.
154	else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
155	altmin=2500.+(40.+ylat)*125.
156	else
157	altmin=2500.
158	endif
159	endif
160
161	do ix=0,nxmin1
162
163	! 1) Calculation of friction velocity
164	!************************************
165
166	ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
167	td2(ix,jy,1,n),surfstr(ix,jy,1,n))
168	if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
169
170	! 2) Calculation of inverse Obukhov length scale
171	!***********************************************
172
173	! NCEP version: find first level above ground
174	llev = 0
175	do i=1,nuvz
176	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
177	end do
178	llev = llev+1
179	if (llev.gt.nuvz) llev = nuvz-1
180	! NCEP version
181
182	! calculate inverse Obukhov length scale with tth(llev)
183	ol=obukhov_gfs(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
184	tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akz(llev))
185	if (ol.ne.0.) then
186	oli(ix,jy,1,n)=1./ol
187	else
188	oli(ix,jy,1,n)=99999.
189	endif
190
191
192	! 3) Calculation of convective velocity scale and mixing height
193	!**************************************************************
194
195	do i=1,nuvz
196	ulev(i)=uuh(ix,jy,i)
197	vlev(i)=vvh(ix,jy,i)
198	ttlev(i)=tth(ix,jy,i,n)
199	qvlev(i)=qvh(ix,jy,i,n)
200	end do
201
202	! NCEP version hmix has been read in in readwind.f, is therefore not calculated here
203	call richardson_gfs(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
204	ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
205	td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus)
206
207	if(lsubgrid.eq.1) then
208	subsceff=min(excessoro(ix,jy),hmixplus)
209	else
210	subsceff=0
211	endif
212	!
213	! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
214	!
215	hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
216	hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
217	hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
218
219	! 4) Calculation of dry deposition velocities
220	!********************************************
221
222	if (DRYDEP) then
223	! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
224	! windspeed
225	z0(7)=0.016ustar(ix,jy,1,n)ustar(ix,jy,1,n)/ga
226
227	! Calculate relative humidity at surface
228	!***************************************
229	rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
230
231	call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
232	tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
233	ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
234	sd(ix,jy,1,n),vd)
235
236	do i=1,nspec
237	vdep(ix,jy,i,n)=vd(i)
238	end do
239
240	endif
241
242	!******************************************************
243	! Calculate height of thermal tropopause (Hoinka, 1997)
244	!******************************************************
245
246	! 1) Calculate altitudes of NCEP model levels
247	!*********************************************
248
249	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
250	ps(ix,jy,1,n))
251	pold=ps(ix,jy,1,n)
252	zold=0.
253	do kz=llev,nuvz
254	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers
255	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
256
257	if (abs(tv-tvold).gt.0.2) then
258	zlev(kz)=zold+constlog(pold/pint)(tv-tvold)/log(tv/tvold)
259	else
260	zlev(kz)=zold+constlog(pold/pint)tv
261	endif
262	tvold=tv
263	pold=pint
264	zold=zlev(kz)
265	end do
266
267	! 2) Define a minimum level kzmin, from which upward the tropopause is
268	! searched for. This is to avoid inversions in the lower troposphere
269	! to be identified as the tropopause
270	!************************************************************************
271
272	do kz=llev,nuvz
273	if (zlev(kz).ge.altmin) then
274	kzmin=kz
275	goto 45
276	endif
277	end do
278	45 continue
279
280	! 3) Search for first stable layer above minimum height that fulfills the
281	! thermal tropopause criterion
282	!************************************************************************
283
284	do kz=kzmin,nuvz
285	do lz=kz+1,nuvz
286	if ((zlev(lz)-zlev(kz)).gt.2000.) then
287	if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
288	(zlev(lz)-zlev(kz))).lt.0.002) then
289	tropopause(ix,jy,1,n)=zlev(kz)
290	goto 51
291	endif
292	goto 50
293	endif
294	end do
295	50 continue
296	end do
297	51 continue
298
299
300	end do
301	end do
302
303
304	! Calculation of potential vorticity on 3-d grid
305	!***********************************************
306
307	call calcpv(n,uuh,vvh,pvh)
308
309
310	end subroutine calcpar_gfs

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

Download in other formats:

Original Format