Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

readwind_timeav.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago
sources for flexwrf v3.1
File size: 16.5 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23
24	subroutine readwind_timeav(indj,n,uuh,vvh,wwh)
25	!**********************************************************************
26	! *
27	! TRAJECTORY MODEL SUBROUTINE READWIND *
28	! *
29	!**********************************************************************
30	! *
31	! April 2012, J. Brioude: This routine handles the difference in time
32	! for time-average fields.
33	!**********************************************************************
34	! *
35	! Note: This is the FLEXPART_WRF version of subroutine readwind. *
36	! The met fields are read from WRF netcdf output files. *
37	! There are many differences from the FLEXPART version. *
38	! *
39	! DESCRIPTION: *
40	! *
41	! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE *
42	! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE *
43	! *
44	! INPUT: *
45	! indj indicates number of the wind field to be read in *
46	! n temporal index for meteorological fields (1 to 3)*
47	! *
48	! IMPORTANT VARIABLES FROM COMMON BLOCK: *
49	! *
50	! wfname File name of data to be read in *
51	! nx,ny,nuvz,nwz expected field dimensions *
52	! nlev_ec number of "T-grid" vertical levels wwf model *
53	! (the unstaggered "bottom_top" dimension) *
54	! uu,vv,ww wind fields *
55	! tt,qv temperature and specific humidity *
56	! ps surface pressure *
57	! *
58	!**********************************************************************
59	!
60
61	! include 'includepar'
62	! include 'includecom'
63	use par_mod
64	use com_mod
65
66	! subr arguments
67	integer :: indj, n
68
69	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
70	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
71	real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
72	real(kind=4) :: uuh2(0:nxmax-1,0:nymax-1,nuvzmax)
73	real(kind=4) :: vvh2(0:nxmax-1,0:nymax-1,nuvzmax)
74	real(kind=4) :: wwh2(0:nxmax-1,0:nymax-1,nwzmax)
75
76	! local variables
77	integer,parameter :: ndims_max=4
78
79	integer :: i, idiagaa, ierr, ifn, itime
80	integer :: iduma,indj2
81	integer :: j, jhhmmss, jyyyymmdd
82	integer :: k, kbgn
83	integer :: lendim(ndims_max), lendim_exp(ndims_max), &
84	lendim_max(ndims_max)
85	integer :: levdiff2,deltat,deltat2
86	integer :: ndims, ndims_exp
87	integer :: n_west_east, n_south_north, n_bottom_top
88	integer :: m_grid_id_dum, m_parent_grid_id_dum, &
89	m_parent_grid_ratio_dum, &
90	i_parent_start_dum, j_parent_start_dum, &
91	map_proj_id_dum, &
92	ext_scalar,pbl_physics,mp_physics_dum
93
94	real :: dx_met, dy_met
95	real :: duma, dumb, dumc, dumd, dume
96	real :: dumdz
97	! real(kind=4) :: dumarray_aa(nwzmax+1)
98	! real(kind=4) :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1)
99	real :: dumarray_aa(nwzmax+1)
100	real :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1)
101	real :: ewater_mb, esatwater_mb
102	real :: ew ! this is an external function
103	real :: map_stdlon_dum, map_truelat1_dum, map_truelat2_dum
104	real :: pint
105	real :: toler
106
107	! real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
108	real :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
109	real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
110
111	real(kind=dp) :: jul,juldate
112
113	character(len=160) :: fnamenc, varname,fnamenc2
114
115	logical :: hflswitch
116
117	!
118	! get grid info from the wrf netcdf file
119	! and check it for consistency against values from gridcheck
120	!
121
122	! print*,'entering timeav'
123
124	fnamenc = path(2)(1:length(2))//wfname(indj)
125	idiagaa = 0
126
127	call read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, &
128	n_west_east, n_south_north, n_bottom_top, &
129	dx_met, dy_met, &
130	m_grid_id_dum, m_parent_grid_id_dum, m_parent_grid_ratio_dum, &
131	i_parent_start_dum, j_parent_start_dum, &
132	map_proj_id_dum, map_stdlon_dum, &
133	map_truelat1_dum, map_truelat2_dum, &
134	ext_scalar,pbl_physics,mp_physics_dum )
135	if (ierr .ne. 0) then
136	write(*,9100) 'error getting gridinfor for met file', fnamenc
137	stop
138	end if
139
140	9100 format( / '*** readwind -- ', a )
141	9110 format( / '*** readwind -- ', a, 1x, i8 / &
142	'file = ', a )
143	9120 format( / '*** readwind -- ', a, 2(1x,i8) / &
144	'file = ', a )
145	9130 format( / '*** readwind -- ', a, 3(1x,i8) / &
146	'file = ', a )
147	9115 format( / '*** readwind -- ', a / a, 1x, i8 / &
148	'file = ', a )
149	9125 format( / '*** readwind -- ', a / a, 2(1x,i8) / &
150	'file = ', a )
151	9135 format( / '*** readwind -- ', a / a, 3(1x,i8) / &
152	'file = ', a )
153
154	toler = 2.0e-7
155
156	if (nx .ne. n_west_east) then
157	write(*,9100) 'nx not consistent', fnamenc
158	stop
159	end if
160	if (ny .ne. n_south_north) then
161	write(*,9100) 'ny not consistent', fnamenc
162	stop
163	end if
164	if (nlev_ec .ne. n_bottom_top) then
165	write(*,9100) 'nlev_ec not consistent', fnamenc
166	stop
167	end if
168	if (nwz .ne. n_bottom_top+1) then
169	write(*,9100) 'nwz not consistent', fnamenc
170	stop
171	end if
172	! if (nuvz .ne. n_bottom_top+1) then
173	! write(*,9100) 'nuvz not consistent', fnamenc
174	! stop
175	! end if
176
177	if (m_grid_id(0) .ne. m_grid_id_dum) then
178	write(*,9100) 'm_grid_id not consistent', fnamenc
179	write(,) m_grid_id(0), m_grid_id_dum
180	stop
181	end if
182	if (m_parent_grid_id(0) .ne. m_parent_grid_id_dum) then
183	write(*,9100) 'm_parent_grid_id not consistent', fnamenc
184	stop
185	end if
186	if (m_parent_grid_ratio(0) .ne. m_parent_grid_ratio_dum) then
187	write(*,9100) 'm_parent_grid_ratio not consistent', fnamenc
188	stop
189	end if
190	if (i_parent_start(0) .ne. i_parent_start_dum) then
191	write(*,9100) 'i_parent_start not consistent', fnamenc
192	stop
193	end if
194	if (j_parent_start(0) .ne. j_parent_start_dum) then
195	write(*,9100) 'j_parent_start not consistent', fnamenc
196	stop
197	end if
198
199	if (abs(dx - dx_met) .gt. toler*abs(dx)) then
200	write(*,9100) 'dx not consistent', fnamenc
201	stop
202	end if
203	if (abs(dy - dy_met) .gt. toler*abs(dy)) then
204	write(*,9100) 'dy not consistent', fnamenc
205	stop
206	end if
207
208
209	! READ THE FIRST FILE
210	if (ldirect.eq.1) indj2=indj
211	if (ldirect.eq.-1) indj2=indj
212	deltat=wfdt(indj2)
213	fnamenc = path(2)(1:length(2))//wfname(indj2)
214	! locate the date/time in the file
215	itime = 0
216	1100 itime = itime + 1
217	call read_ncwrfout_1datetime( ierr, fnamenc, &
218	itime, jyyyymmdd, jhhmmss )
219	if (ierr .eq. -1) then
220	write(*,9100) 'error reading time from met file', fnamenc
221	stop
222	else if (ierr .ne. 0) then
223	write(*,9125) 'unable to locate date/time in met file', &
224	'indj, itime =', indj2, itime, fnamenc
225	stop
226	else
227	jul = juldate( jyyyymmdd, jhhmmss )
228	duma = (jul-bdate)*86400.
229	iduma = nint(duma)
230	if (iduma .ne. wftime(indj2)) goto 1100
231	end if
232	if (option_verbose.eq.1) then
233	write(,)
234	write(,) 'readwind processing wrfout file ='
235	write(,) fnamenc
236	write(,) 'itime, ymd, hms =', itime, jyyyymmdd, jhhmmss
237	endif
238	kbgn = 1 + add_sfc_level
239
240	if (wind_option.eq.1) varname = 'AVGFLX_RUM'
241	do i = 1, ndims_max
242	lendim_exp(i) = 0
243	lendim_max(i) = 1
244	end do
245	lendim_exp(1) = nx+1
246	lendim_max(1) = nxmax
247	lendim_exp(2) = ny
248	lendim_max(2) = nymax
249	lendim_exp(3) = nuvz-add_sfc_level
250	lendim_max(3) = nuvzmax
251	ndims_exp = 4
252	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
253	varname, uuh(0,0,kbgn), &
254	itime, &
255	ndims, ndims_exp, ndims_max, &
256	lendim, lendim_exp, lendim_max )
257	if (ierr .ne. 0) then
258	write(*,9100) 'error doing ncread of U', fnamenc
259	if (wind_option.le.0) print*,'you asked snapshot winds'
260	if (wind_option.eq.1) print*,'you asked mean winds'
261	print*,'change wind_option'
262
263	stop
264	end if
265
266
267	! v wind velocity
268	! the wrf output file contains (nuvz-add_sfc_level) levels
269	! read the data into k=kbgn,nuvz
270	! (interpolate it from "V-grid" to "T-grid" later)
271	if (wind_option.eq.1) varname = 'AVGFLX_RVM'
272	lendim_exp(1) = nx
273	lendim_max(1) = nxmax
274	lendim_exp(2) = ny+1
275	lendim_max(2) = nymax
276	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
277	varname, vvh(0,0,kbgn), &
278	itime, &
279	ndims, ndims_exp, ndims_max, &
280	lendim, lendim_exp, lendim_max )
281	if (ierr .ne. 0) then
282	write(*,9100) 'error doing ncread of V', fnamenc
283	if (wind_option.eq.0) print*,'you asked snapshot winds'
284	if (wind_option.eq.1) print*,'you asked mean winds'
285	print*,'change wind_option'
286	stop
287	end if
288
289
290	! w wind velocity
291	! this is on the "W-grid", and
292	! the wrf output file contains nwz levels, so no shifting needed
293	if (wind_option.eq.1) varname = 'AVGFLX_WWM'
294	! print*,'varname',varname
295	lendim_exp(1) = nx
296	lendim_max(1) = nxmax
297	lendim_exp(2) = ny
298	lendim_max(2) = nymax
299	lendim_exp(3) = nwz
300	lendim_max(3) = nwzmax
301	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
302	varname, wwh, &
303	itime, &
304	ndims, ndims_exp, ndims_max, &
305	lendim, lendim_exp, lendim_max )
306	if (ierr .ne. 0) then
307	write(*,9100) 'error doing ncread of W', fnamenc
308	if (wind_option.eq.0) print*,'you asked snapshot winds'
309	if (wind_option.eq.1) print*,'you asked mean winds'
310	print*,'change wind_option'
311	stop
312	end if
313
314	! READ THE SECOND FILE and average
315	! deltat must be equal to deltat2, otherwise the time-average wind cannot be
316	! fixed.
317	! deltat is assumed to be the same than the WRF output.
318
319	if (ldirect.eq.1) indj2=indj+1
320	if (ldirect.eq.-1) indj2=indj-1
321	deltat2=wfdt(indj2)
322	fnamenc = path(2)(1:length(2))//wfname(indj2)
323
324	! print*,'test timeav',deltat,deltat2,indj2,numbwf
325	if (deltat.eq.deltat2 .and. indj2.ge.1 .and.indj2.le.numbwf) then
326	! print*,'reading second file in timeav',deltat,deltat2,indj2
327	! locate the date/time in the file
328	itime = 0
329	1101 itime = itime + 1
330	call read_ncwrfout_1datetime( ierr, fnamenc, &
331	itime, jyyyymmdd, jhhmmss )
332	if (ierr .eq. -1) then
333	write(*,9100) 'error reading time from met file', fnamenc
334	stop
335	else if (ierr .ne. 0) then
336	write(*,9125) 'unable to locate date/time in met file', &
337	'indj, itime =', indj2, itime, fnamenc
338	stop
339	else
340	jul = juldate( jyyyymmdd, jhhmmss )
341	duma = (jul-bdate)*86400.
342	iduma = nint(duma)
343	if (iduma .ne. wftime(indj2)) goto 1101
344	end if
345	if (option_verbose.eq.1) then
346	write(,)
347	write(,) 'readwind processing wrfout file ='
348	write(,) fnamenc
349	write(,) 'itime, ymd, hms =', itime, jyyyymmdd, jhhmmss
350	endif
351	kbgn = 1 + add_sfc_level
352
353	if (wind_option.eq.1) varname = 'AVGFLX_RUM'
354	! print*,'varname',varname
355	do i = 1, ndims_max
356	lendim_exp(i) = 0
357	lendim_max(i) = 1
358	end do
359	lendim_exp(1) = nx+1
360	lendim_max(1) = nxmax
361	lendim_exp(2) = ny
362	lendim_max(2) = nymax
363	lendim_exp(3) = nuvz-add_sfc_level
364	lendim_max(3) = nuvzmax
365	ndims_exp = 4
366	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
367	varname, uuh2(0,0,kbgn), &
368	itime, &
369	ndims, ndims_exp, ndims_max, &
370	lendim, lendim_exp, lendim_max )
371	if (ierr .ne. 0) then
372	write(*,9100) 'error doing ncread of U', fnamenc
373	if (wind_option.le.0) print*,'you asked snapshot winds'
374	if (wind_option.eq.1) print*,'you asked mean winds'
375	print*,'change wind_option'
376
377	stop
378	end if
379
380
381	! v wind velocity
382	! the wrf output file contains (nuvz-add_sfc_level) levels
383	! read the data into k=kbgn,nuvz
384	! (interpolate it from "V-grid" to "T-grid" later)
385	if (wind_option.eq.1) varname = 'AVGFLX_RVM'
386	! print*,'varname',varname
387	lendim_exp(1) = nx
388	lendim_max(1) = nxmax
389	lendim_exp(2) = ny+1
390	lendim_max(2) = nymax
391	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
392	varname, vvh2(0,0,kbgn), &
393	itime, &
394	ndims, ndims_exp, ndims_max, &
395	lendim, lendim_exp, lendim_max )
396	if (ierr .ne. 0) then
397	write(*,9100) 'error doing ncread of V', fnamenc
398	if (wind_option.eq.0) print*,'you asked snapshot winds'
399	if (wind_option.eq.1) print*,'you asked mean winds'
400	print*,'change wind_option'
401	stop
402	end if
403
404
405	! w wind velocity
406	! this is on the "W-grid", and
407	! the wrf output file contains nwz levels, so no shifting needed
408	if (wind_option.eq.1) varname = 'AVGFLX_WWM'
409	! print*,'varname',varname
410	lendim_exp(1) = nx
411	lendim_max(1) = nxmax
412	lendim_exp(2) = ny
413	lendim_max(2) = nymax
414	lendim_exp(3) = nwz
415	lendim_max(3) = nwzmax
416	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
417	varname, wwh2, &
418	itime, &
419	ndims, ndims_exp, ndims_max, &
420	lendim, lendim_exp, lendim_max )
421	if (ierr .ne. 0) then
422	write(*,9100) 'error doing ncread of W', fnamenc
423	if (wind_option.eq.0) print*,'you asked snapshot winds'
424	if (wind_option.eq.1) print*,'you asked mean winds'
425	print*,'change wind_option'
426	stop
427	end if
428
429	! print*,'start modif timeav'
430	do k=1,nuvzmax
431	do j=0,nymax-1
432	do i=0,nxmax-1
433	uuh(i,j,k)=0.5*(uuh(i,j,k)+uuh2(i,j,k))
434	vvh(i,j,k)=0.5*(vvh(i,j,k)+vvh2(i,j,k))
435	enddo
436	enddo
437	enddo
438	do k=1,nwzmax
439	do j=0,nymax-1
440	do i=0,nxmax-1
441	wwh(i,j,k)=0.5*(wwh(i,j,k)+wwh2(i,j,k))
442	enddo
443	enddo
444	enddo
445	! print*,'end modif timeav'
446
447
448	endif ! test on deltat
449
450	! print*,'out of entering timeav'
451
452	return
453	end subroutine readwind_timeav
454

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

Download in other formats:

Original Format