Context Navigation

calcpar.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug

Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago
move license from headers to a different file
Property mode set to `100644`
File size: 9.9 KB

Line
1	subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
2	! i i i o
3	!*****************************************************************************
4	! *
5	! Computation of several boundary layer parameters needed for the *
6	! dispersion calculation and calculation of dry deposition velocities. *
7	! All parameters are calculated over the entire grid. *
8	! *
9	! Author: A. Stohl *
10	! *
11	! 21 May 1995 *
12	! *
13	! ------------------------------------------------------------------ *
14	! Petra Seibert, Feb 2000: *
15	! convection scheme: *
16	! new variables in call to richardson *
17	! *
18	! CHANGE 17/11/2005 Caroline Forster NCEP GFS version *
19	! *
20	! Changes, Bernd C. Krueger, Feb. 2001: *
21	! Variables tth and qvh (on eta coordinates) in common block *
22	! *
23	! Unified ECMWF and GFS builds *
24	! Marian Harustak, 12.5.2017 *
25	! - Merged calcpar and calcpar_gfs into one routine using if-then *
26	! for meteo-type dependent code *
27	!*****************************************************************************
28
29	!*****************************************************************************
30	! *
31	! Variables: *
32	! n temporal index for meteorological fields (1 to 3) *
33	! uuh *
34	! vvh *
35	! pvh *
36	! metdata_format format of metdata (ecmwf/gfs) *
37	! *
38	! Constants: *
39	! *
40	! *
41	! Functions: *
42	! scalev computation of ustar *
43	! obukhov computatio of Obukhov length *
44	! *
45	!*****************************************************************************
46
47	use par_mod
48	use com_mod
49	use class_gribfile
50
51	implicit none
52
53	integer :: metdata_format
54	integer :: n,ix,jy,i,kz,lz,kzmin,llev,loop_start
55	real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
56	real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
57	real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy,akzdummy
58	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
59	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
60	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
61	real,parameter :: const=r_air/ga
62
63	!write(,) 'in calcpar writting snowheight'
64	!***********************************
65	! for test: write out snow depths
66
67	! open(4,file='slandusetest',form='formatted')
68	! do 5 ix=0,nxmin1
69	!5 write (4,*) (sd(ix,jy,1,n),jy=0,nymin1)
70	! close(4)
71
72
73	! Loop over entire grid
74	!**********************
75
76	do jy=0,nymin1
77
78	! Set minimum height for tropopause
79	!**********************************
80
81	ylat=ylat0+real(jy)*dy
82	if ((ylat.ge.-20.).and.(ylat.le.20.)) then
83	altmin = 5000.
84	else
85	if ((ylat.gt.20.).and.(ylat.lt.40.)) then
86	altmin=2500.+(40.-ylat)*125.
87	else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
88	altmin=2500.+(40.+ylat)*125.
89	else
90	altmin=2500.
91	endif
92	endif
93
94	do ix=0,nxmin1
95
96	! 1) Calculation of friction velocity
97	!************************************
98
99	ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
100	td2(ix,jy,1,n),surfstr(ix,jy,1,n))
101	if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
102
103	! 2) Calculation of inverse Obukhov length scale
104	!***********************************************
105
106	if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
107	! NCEP version: find first level above ground
108	llev = 0
109	do i=1,nuvz
110	if (ps(ix,jy,1,n).lt.akz(i)) llev=i
111	end do
112	llev = llev+1
113	if (llev.gt.nuvz) llev = nuvz-1
114	! NCEP version
115
116	! calculate inverse Obukhov length scale with tth(llev)
117	ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
118	tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
119	else
120	llev=0
121	ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
122	tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akzdummy,metdata_format)
123	end if
124
125	if (ol.ne.0.) then
126	oli(ix,jy,1,n)=1./ol
127	else
128	oli(ix,jy,1,n)=99999.
129	endif
130
131
132	! 3) Calculation of convective velocity scale and mixing height
133	!**************************************************************
134
135	do i=1,nuvz
136	ulev(i)=uuh(ix,jy,i)
137	vlev(i)=vvh(ix,jy,i)
138	ttlev(i)=tth(ix,jy,i,n)
139	qvlev(i)=qvh(ix,jy,i,n)
140	end do
141
142	if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
143	! NCEP version hmix has been read in in readwind.f, is therefore not calculated here
144	call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
145	ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
146	td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus,metdata_format)
147	else
148	call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
149	ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
150	td2(ix,jy,1,n),hmix(ix,jy,1,n),wstar(ix,jy,1,n),hmixplus,metdata_format)
151	end if
152
153	if(lsubgrid.eq.1) then
154	subsceff=min(excessoro(ix,jy),hmixplus)
155	else
156	subsceff=0.0
157	endif
158	!
159	! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
160	!
161	hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
162	hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
163	hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
164
165	! 4) Calculation of dry deposition velocities
166	!********************************************
167
168	if (DRYDEP) then
169	! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
170	! windspeed
171	z0(7)=0.016ustar(ix,jy,1,n)ustar(ix,jy,1,n)/ga
172
173	! Calculate relative humidity at surface
174	!***************************************
175	rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
176
177	call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
178	tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
179	ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
180	sd(ix,jy,1,n),vd)
181
182	do i=1,nspec
183	vdep(ix,jy,i,n)=vd(i)
184	end do
185
186	endif
187
188	!******************************************************
189	! Calculate height of thermal tropopause (Hoinka, 1997)
190	!******************************************************
191
192	! 1) Calculate altitudes of model levels
193	!***************************************
194
195	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
196	ps(ix,jy,1,n))
197	pold=ps(ix,jy,1,n)
198	zold=0.
199	if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
200	loop_start=2
201	else
202	loop_start=llev
203	end if
204	do kz=loop_start,nuvz
205	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) ! pressure on model layers
206	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
207
208	if (abs(tv-tvold).gt.0.2) then
209	zlev(kz)=zold+constlog(pold/pint)(tv-tvold)/log(tv/tvold)
210	else
211	zlev(kz)=zold+constlog(pold/pint)tv
212	endif
213	tvold=tv
214	pold=pint
215	zold=zlev(kz)
216	end do
217
218	! 2) Define a minimum level kzmin, from which upward the tropopause is
219	! searched for. This is to avoid inversions in the lower troposphere
220	! to be identified as the tropopause
221	!************************************************************************
222
223	if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
224	loop_start=1
225	else
226	loop_start=llev
227	end if
228
229	do kz=loop_start,nuvz
230	if (zlev(kz).ge.altmin) then
231	kzmin=kz
232	goto 45
233	endif
234	end do
235	45 continue
236
237	! 3) Search for first stable layer above minimum height that fulfills the
238	! thermal tropopause criterion
239	!************************************************************************
240
241	do kz=kzmin,nuvz
242	do lz=kz+1,nuvz
243	if ((zlev(lz)-zlev(kz)).gt.2000.) then
244	if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
245	(zlev(lz)-zlev(kz))).lt.0.002) then
246	tropopause(ix,jy,1,n)=zlev(kz)
247	goto 51
248	endif
249	goto 50
250	endif
251	end do
252	50 continue
253	end do
254	51 continue
255
256
257	end do
258	end do
259
260	! Calculation of potential vorticity on 3-d grid
261	!***********************************************
262
263	call calcpv(n,uuh,vvh,pvh)
264
265
266	end subroutine calcpar

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

Download in other formats:

Original Format