Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

readwind.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 47.9 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(indj,n,uuh,vvh,wwh,divh)
25	!**********************************************************************
26	! *
27	! TRAJECTORY MODEL SUBROUTINE READWIND *
28	! *
29	!**********************************************************************
30	! *
31	! AUTHOR: G. WOTAWA *
32	! DATE: 1997-08-05 *
33	! LAST UPDATE: 2000-10-17, Andreas Stohl *
34	! *
35	! Bernd C. Krueger, Feb. 2001: Variables tth and qvh *
36	! (on eta coordinates) in common block *
37	! *
38	! Oct-Dec, 2005: R. Easter. Major changes for WRF. *
39	! 06-nov-2005 rce - change uuh,vvh dimension back to original *
40	! 16-nov-2005 rce - zzh is shifted like pph,tth *
41	! *
42	! 11-June-2007, W.WANG -- read TKE, change ndims_exp=1 for P_TOP
43	! 19-Oct -2007, read tcc, RAINC, RAINNC, CLDFRA,W0AVG
44	! Note RAINC, RAINNC are accumulated prec
45	! Dec 2011, J Brioude: modifications, notably for the mean wind
46	!
47	! D. Arnold May 2012: quick fix - for CLDFRA
48	! CLDFRA was for sub-grid variability of precipitation
49	! within a cell from Hertel et al., 1995 (grid cells of approx. 150 km)
50	! WRF will be typically be run at higher resolutions, therefore
51	! we simply skip this condition and consider the maximum fraction
52	! no more reading of CLDFRA needed, initialized to 0.
53	! Same modifications for the nested domains!
54
55	!**********************************************************************
56	! *
57	! Note: This is the FLEXPART_WRF version of subroutine readwind. *
58	! The met fields are read from WRF netcdf output files. *
59	! There are many differences from the FLEXPART version. *
60	! *
61	! DESCRIPTION: *
62	! *
63	! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE *
64	! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE *
65	! *
66	! INPUT: *
67	! indj indicates number of the wind field to be read in *
68	! n temporal index for meteorological fields (1 to 3)*
69	! *
70	! IMPORTANT VARIABLES FROM COMMON BLOCK: *
71	! *
72	! wfname File name of data to be read in *
73	! nx,ny,nuvz,nwz expected field dimensions *
74	! nlev_ec number of "T-grid" vertical levels wwf model *
75	! (the unstaggered "bottom_top" dimension) *
76	! uu,vv,ww wind fields *
77	! tt,qv temperature and specific humidity *
78	! ps surface pressure *
79	! *
80	!**********************************************************************
81	!
82
83	! include 'includepar'
84	! include 'includecom'
85	use par_mod
86	use com_mod
87
88	! subr arguments
89	integer :: indj, n
90
91	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
92	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
93	real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
94	! real(kind=4) :: urot(0:nxmax-1,0:nymax-1,nuvzmax)
95	! real(kind=4) :: vrot(0:nxmax-1,0:nymax-1,nuvzmax)
96	real(kind=4) :: divh(0:nxmax-1,0:nymax-1,nuvzmax)
97	real(kind=4) :: mu(0:nxmax-1,0:nymax-1,1),mu2
98	real(kind=4) :: mub(0:nxmax-1,0:nymax-1,1)
99	! real :: utrue1,vtrue1,utrue2,vtrue2,dumy
100	! real(kind=4) :: m_u(0:nxmax-1,0:nymax-1,1)
101	! real(kind=4) :: m_v(0:nxmax-1,0:nymax-1,1)
102	! real(kind=4) :: m_w(0:nxmax-1,0:nymax-1,1)
103
104	! real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
105	! real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
106	! real :: urot(0:nxmax-1,0:nymax-1,nuvzmax)
107	! real :: vrot(0:nxmax-1,0:nymax-1,nuvzmax)
108	! real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
109	! real :: divh(0:nxmax-1,0:nymax-1,nuvzmax)
110	! real :: mu(0:nxmax-1,0:nymax-1,1),mu2
111	! real :: mub(0:nxmax-1,0:nymax-1,1)
112	! real :: m_u(0:nxmax-1,0:nymax-1,1)
113	! real :: m_v(0:nxmax-1,0:nymax-1,1)
114	! real :: m_x(0:nxmax-1,0:nymax-1,1)
115	! real :: m_y(0:nxmax-1,0:nymax-1,1)
116	! real :: m_z(0:nxmax-1,0:nymax-1,1)
117
118	! local variables
119	integer,parameter :: ndims_max=4
120
121	integer :: i, idiagaa, ierr, ifn, itime
122	integer :: iduma
123	integer :: j, jhhmmss, jyyyymmdd
124	integer :: k, kbgn
125	integer :: lendim(ndims_max), lendim_exp(ndims_max), &
126	lendim_max(ndims_max)
127	integer :: levdiff2
128	integer :: ndims, ndims_exp
129	integer :: n_west_east, n_south_north, n_bottom_top
130	integer :: m_grid_id_dum, m_parent_grid_id_dum, &
131	m_parent_grid_ratio_dum, &
132	i_parent_start_dum, j_parent_start_dum, &
133	map_proj_id_dum, &
134	ext_scalar,pbl_physics,mp_physics_dum
135
136	real :: dx_met, dy_met
137	real :: duma, dumb, dumc, dumd, dume
138	real :: dumdz
139	! real(kind=4) :: dumarray_aa(nwzmax+1)
140	! real(kind=4) :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1)
141	real :: dumarray_aa(nwzmax+1)
142	real :: dumarray_pp(0:nxmax-1,0:nymax-1,nwzmax+1)
143	real :: ewater_mb, esatwater_mb
144	real :: ew ! this is an external function
145	real :: map_stdlon_dum, map_truelat1_dum, map_truelat2_dum
146	real :: pint
147	real :: toler
148
149	! real(kind=4) :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
150	real :: ewss(0:nxmax-1,0:nymax-1),nsss(0:nxmax-1,0:nymax-1)
151	real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
152
153	real(kind=dp) :: jul,juldate
154
155	character(len=160) :: fnamenc, varname,fnamenc2
156
157	logical :: hflswitch
158
159	!
160	! get grid info from the wrf netcdf file
161	! and check it for consistency against values from gridcheck
162	!
163	fnamenc = path(2)(1:length(2))//wfname(indj)
164	idiagaa = 0
165
166	call read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, &
167	n_west_east, n_south_north, n_bottom_top, &
168	dx_met, dy_met, &
169	m_grid_id_dum, m_parent_grid_id_dum, m_parent_grid_ratio_dum, &
170	i_parent_start_dum, j_parent_start_dum, &
171	map_proj_id_dum, map_stdlon_dum, &
172	map_truelat1_dum, map_truelat2_dum, &
173	ext_scalar,pbl_physics,mp_physics_dum )
174	if (ierr .ne. 0) then
175	write(*,9100) 'error getting gridinfor for met file', fnamenc
176	stop
177	end if
178
179	9100 format( / '*** readwind -- ', a )
180	9110 format( / '*** readwind -- ', a, 1x, i8 / &
181	'file = ', a )
182	9120 format( / '*** readwind -- ', a, 2(1x,i8) / &
183	'file = ', a )
184	9130 format( / '*** readwind -- ', a, 3(1x,i8) / &
185	'file = ', a )
186	9115 format( / '*** readwind -- ', a / a, 1x, i8 / &
187	'file = ', a )
188	9125 format( / '*** readwind -- ', a / a, 2(1x,i8) / &
189	'file = ', a )
190	9135 format( / '*** readwind -- ', a / a, 3(1x,i8) / &
191	'file = ', a )
192
193	toler = 2.0e-7
194
195	if (nx .ne. n_west_east) then
196	write(*,9100) 'nx not consistent', fnamenc
197	stop
198	end if
199	if (ny .ne. n_south_north) then
200	write(*,9100) 'ny not consistent', fnamenc
201	stop
202	end if
203	if (nlev_ec .ne. n_bottom_top) then
204	write(*,9100) 'nlev_ec not consistent', fnamenc
205	stop
206	end if
207	if (nwz .ne. n_bottom_top+1) then
208	write(*,9100) 'nwz not consistent', fnamenc
209	stop
210	end if
211	! if (nuvz .ne. n_bottom_top+1) then
212	! write(*,9100) 'nuvz not consistent', fnamenc
213	! stop
214	! end if
215
216	if (m_grid_id(0) .ne. m_grid_id_dum) then
217	write(*,9100) 'm_grid_id not consistent', fnamenc
218	write(,) m_grid_id(0), m_grid_id_dum
219	stop
220	end if
221	if (m_parent_grid_id(0) .ne. m_parent_grid_id_dum) then
222	write(*,9100) 'm_parent_grid_id not consistent', fnamenc
223	stop
224	end if
225	if (m_parent_grid_ratio(0) .ne. m_parent_grid_ratio_dum) then
226	write(*,9100) 'm_parent_grid_ratio not consistent', fnamenc
227	stop
228	end if
229	if (i_parent_start(0) .ne. i_parent_start_dum) then
230	write(*,9100) 'i_parent_start not consistent', fnamenc
231	stop
232	end if
233	if (j_parent_start(0) .ne. j_parent_start_dum) then
234	write(*,9100) 'j_parent_start not consistent', fnamenc
235	stop
236	end if
237
238	if (abs(dx - dx_met) .gt. toler*abs(dx)) then
239	write(*,9100) 'dx not consistent', fnamenc
240	stop
241	end if
242	if (abs(dy - dy_met) .gt. toler*abs(dy)) then
243	write(*,9100) 'dy not consistent', fnamenc
244	stop
245	end if
246
247	! locate the date/time in the file
248	itime = 0
249	1100 itime = itime + 1
250	call read_ncwrfout_1datetime( ierr, fnamenc, &
251	itime, jyyyymmdd, jhhmmss )
252	if (ierr .eq. -1) then
253	write(*,9100) 'error reading time from met file', fnamenc
254	stop
255	else if (ierr .ne. 0) then
256	write(*,9125) 'unable to locate date/time in met file', &
257	'indj, itime =', indj, itime, fnamenc
258	stop
259	else
260	jul = juldate( jyyyymmdd, jhhmmss )
261	duma = (jul-bdate)*86400.
262	iduma = nint(duma)
263	if (iduma .ne. wftime(indj)) goto 1100
264	end if
265	if (option_verbose.eq.1) then
266	write(,)
267	write(,) 'readwind processing wrfout file ='
268	write(,) fnamenc
269	write(,) 'itime, ymd, hms =', itime, jyyyymmdd, jhhmmss
270
271	endif
272	! read eta_w_wrf, eta_u_wrf, p_top_wrf, ylat2d, xlon2d from the
273	! netcdf wrfout file and compare to those from the 1st met. file
274
275	varname = 'ZNW'
276	do i = 1, ndims_max
277	lendim_exp(i) = 0
278	lendim_max(i) = 1
279	end do
280	lendim_exp(1) = nwz
281	lendim_max(1) = nwzmax+1
282	ndims_exp = 2
283	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
284	varname, dumarray_aa, &
285	itime, &
286	ndims, ndims_exp, ndims_max, &
287	lendim, lendim_exp, lendim_max )
288	if (ierr .ne. 0) then
289	fnamenc2='wrfout_d03_zn.nc'
290	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, &
291	varname, dumarray_aa, &
292	itime, &
293	ndims, ndims_exp, ndims_max, &
294	lendim, lendim_exp, lendim_max )
295
296	if (ierr .ne. 0) then
297	write(*,9100) 'error doing ncread of ZNW', fnamenc
298	stop
299	end if
300	end if
301	do k = 1, nwz
302	if (abs(eta_w_wrf(k) - dumarray_aa(k)) &
303	.gt. toler*abs(eta_w_wrf(k))) then
304	write(*,9100) 'eta_w_wrf not consistent', fnamenc
305	stop
306	end if
307	end do
308
309	varname = 'ZNU'
310	lendim_exp(1) = nwz-1
311	lendim_max(1) = nwzmax+1
312	ndims_exp = 2
313	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
314	varname, dumarray_aa, &
315	itime, &
316	ndims, ndims_exp, ndims_max, &
317	lendim, lendim_exp, lendim_max )
318	if (ierr .ne. 0) then
319	fnamenc2='wrfout_d03_zn.nc'
320	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, &
321	varname, dumarray_aa, &
322	itime, &
323	ndims, ndims_exp, ndims_max, &
324	lendim, lendim_exp, lendim_max )
325	if (ierr .ne. 0) then
326
327	write(*,9100) 'error doing ncread of ZNU', fnamenc
328	stop
329	end if
330	end if
331	do k = 1, nwz-1
332	if (abs(eta_u_wrf(k) - dumarray_aa(k)) &
333	.gt. toler*abs(eta_u_wrf(k))) then
334	write(*,9100) 'eta_u_wrf not consistent', fnamenc
335	stop
336	end if
337	end do
338
339	! varname = 'P_TOP'
340	! lendim_exp(1) = 1
341	! lendim_max(1) = 6
342	! ndims_exp = 2
343	! if (ext_scalar .lt. 0) ndims_exp = 1
344	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
345	! varname, duma, &
346	! itime, &
347	! ndims, ndims_exp, ndims_max, &
348	! lendim, lendim_exp, lendim_max )
349	! if (ierr .ne. 0) then
350	! write(*,9100) 'error doing ncread of P_TOP', fnamenc
351	! stop
352	! end if
353	! if (abs(p_top_wrf - duma) .gt. toler*abs(p_top_wrf)) then
354	! write(*,9100) 'p_top_wrf not consistent', fnamenc
355	! stop
356	! end if
357
358	varname = 'XLAT'
359	lendim_exp(1) = nx
360	lendim_max(1) = nxmax
361	lendim_exp(2) = ny
362	lendim_max(2) = nymax
363	ndims_exp = 3
364	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
365	varname, dumarray_pp, &
366	itime, &
367	ndims, ndims_exp, ndims_max, &
368	lendim, lendim_exp, lendim_max )
369	if (ierr .ne. 0) then
370	write(,)
371	write(*,9100) 'error doing ncread of XLAT', fnamenc
372	stop
373	end if
374	toler = 1.0e-6
375	do j = 0, ny-1
376	do i = 0, nx-1
377	if (abs(ylat2d(i,j) - dumarray_pp(i,j,1)) .gt. &
378	toler*abs(ylat2d(i,j))) then
379	write(*,9100) 'ylat2d not consistent', fnamenc
380	write(*,'(a,2i5,2f16.6)') 'i,j,ylats =', i, j, &
381	ylat2d(i,j), dumarray_pp(i,j,1)
382	stop
383	end if
384	end do
385	end do
386
387	varname = 'XLONG'
388	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
389	varname, dumarray_pp, &
390	itime, &
391	ndims, ndims_exp, ndims_max, &
392	lendim, lendim_exp, lendim_max )
393	if (ierr .ne. 0) then
394	write(,)
395	write(*,9100) 'error doing ncread of XLONG', fnamenc
396	stop
397	end if
398	do j = 0, ny-1
399	do i = 0, nx-1
400	if (abs(xlon2d(i,j) - dumarray_pp(i,j,1)) .gt. &
401	toler*abs(xlon2d(i,j))) then
402	write(*,9100) 'xlon2d not consistent', fnamenc
403	write(*,'(a,2i5,2f16.6)') 'i,j,xlons =', i, j, &
404	xlon2d(i,j), dumarray_pp(i,j,1)
405	stop
406	end if
407	end do
408	end do
409
410
411	!
412	!
413	! now read the data fields for current time
414	! the following are read from ecmwf met files
415	! U VELOCITY
416	! V VELOCITY
417	! W VELOCITY
418	! TEMPERATURE
419	! SPEC. HUMIDITY
420	! SURF. PRESS.
421	! SEA LEVEL PRESS.
422	! 10 M U VELOCITY
423	! 10 M V VELOCITY
424	! 2 M TEMPERATURE
425	! 2 M DEW POINT
426	! SNOW DEPTH
427	! CLOUD COVER
428	! LARGE SCALE PREC.
429	! CONVECTIVE PREC.
430	! SENS. HEAT FLUX
431	! SOLAR RADIATION
432	! EW SURFACE STRESS
433	! NS SURFACE STRESS
434	! ECMWF OROGRAPHY
435	! STANDARD DEVIATION OF OROGRAPHY
436	! ECMWF LAND SEA MASK
437	!
438	!
439	hflswitch=.false.
440	strswitch=.false.
441	levdiff2=nlev_ec-nwz+1
442
443	kbgn = 1 + add_sfc_level
444	! at this point
445	! if add_sfc_level=1, then nuvz=nwz and kbgn=2
446	! if add_sfc_level=0, then nuvz=nwz-1 and kbgn=1
447
448	! u wind velocity
449	! the wrf output file contains (nuvz-add_sfc_level) levels
450	! read the data into k=kbgn,nuvz
451	! (interpolate it from "U-grid" to "T-grid" later)
452	if (wind_option.le.0) varname = 'U'
453	if (wind_option.eq.1) varname = 'AVGFLX_RUM'
454	if (wind_option.eq.2) varname = 'U'
455
456	do i = 1, ndims_max
457	lendim_exp(i) = 0
458	lendim_max(i) = 1
459	end do
460	lendim_exp(1) = nx+1
461	lendim_max(1) = nxmax
462	lendim_exp(2) = ny
463	lendim_max(2) = nymax
464	lendim_exp(3) = nuvz-add_sfc_level
465	lendim_max(3) = nuvzmax
466	ndims_exp = 4
467	if (time_option.eq.0) then
468	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
469	varname, uuh(0,0,kbgn), &
470	itime, &
471	ndims, ndims_exp, ndims_max, &
472	lendim, lendim_exp, lendim_max )
473	if (ierr .ne. 0) then
474	write(*,9100) 'error doing ncread of U', fnamenc
475	if (wind_option.le.0) print*,'you asked snapshot winds'
476	if (wind_option.eq.1) print*,'you asked mean winds'
477	print*,'change wind_option'
478
479	stop
480	end if
481
482	endif !test on time_option
483
484	! v wind velocity
485	! the wrf output file contains (nuvz-add_sfc_level) levels
486	! read the data into k=kbgn,nuvz
487	! (interpolate it from "V-grid" to "T-grid" later)
488	if (wind_option.le.0) varname = 'V'
489	if (wind_option.eq.1) varname = 'AVGFLX_RVM'
490	if (wind_option.eq.2) varname = 'V'
491	lendim_exp(1) = nx
492	lendim_max(1) = nxmax
493	lendim_exp(2) = ny+1
494	lendim_max(2) = nymax
495	if (time_option.eq.0) then
496	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
497	varname, vvh(0,0,kbgn), &
498	itime, &
499	ndims, ndims_exp, ndims_max, &
500	lendim, lendim_exp, lendim_max )
501	if (ierr .ne. 0) then
502	write(*,9100) 'error doing ncread of V', fnamenc
503	if (wind_option.eq.0) print*,'you asked snapshot winds'
504	if (wind_option.eq.1) print*,'you asked mean winds'
505	print*,'change wind_option'
506	stop
507	end if
508
509	endif !test on time_option
510
511	! w wind velocity
512	! this is on the "W-grid", and
513	! the wrf output file contains nwz levels, so no shifting needed
514	if (wind_option.le.0) varname = 'W'
515	if (wind_option.eq.1) varname = 'AVGFLX_WWM'
516	if (wind_option.eq.2) varname = 'WW'
517	! print*,'varname',varname
518	lendim_exp(1) = nx
519	lendim_max(1) = nxmax
520	lendim_exp(2) = ny
521	lendim_max(2) = nymax
522	lendim_exp(3) = nwz
523	lendim_max(3) = nwzmax
524	if (time_option.eq.0) then
525	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
526	varname, wwh, &
527	itime, &
528	ndims, ndims_exp, ndims_max, &
529	lendim, lendim_exp, lendim_max )
530	if (ierr .ne. 0) then
531	write(*,9100) 'error doing ncread of W', fnamenc
532	if (wind_option.eq.0) print*,'you asked snapshot winds'
533	if (wind_option.eq.1) print*,'you asked mean winds'
534	print*,'change wind_option'
535	stop
536	end if
537
538	endif !test on time_option
539
540	! pressure - read base state and perturbation pressure,
541	! then combine
542	! the wrf output file contains (nuvz-add_sfc_level) levels
543	! read the data into k=kbgn,nuvz
544	varname = 'PB'
545	lendim_exp(3) = nuvz-add_sfc_level
546	lendim_max(3) = nuvzmax
547	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
548	varname, pph(0,0,kbgn,n), &
549	itime, &
550	ndims, ndims_exp, ndims_max, &
551	lendim, lendim_exp, lendim_max )
552	if (ierr .ne. 0) then
553	write(*,9100) 'error doing ncread of PB', fnamenc
554	stop
555	end if
556
557	varname = 'P'
558	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
559	varname, dumarray_pp(0,0,kbgn), &
560	itime, &
561	ndims, ndims_exp, ndims_max, &
562	lendim, lendim_exp, lendim_max )
563	if (ierr .ne. 0) then
564	write(*,9100) 'error doing ncread of P', fnamenc
565	stop
566	end if
567
568	do k = kbgn, nuvz
569	do j = 0, nymin1
570	do i = 0, nxmin1
571	pph(i,j,k,n) = pph(i,j,k,n) + dumarray_pp(i,j,k)
572	end do
573	end do
574	end do
575
576
577	! height - read base state and perturbation geopotential,
578	! then combine and divide by gravity
579	! these are on the "W-grid", and
580	! the wrf output file contains nwz levels
581	! shift them also so they will be consistent with pph
582	varname = 'PHB'
583	lendim_exp(3) = nwz
584	lendim_max(3) = nwzmax+1
585	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
586	varname, zzh(0,0,kbgn,n), &
587	itime, &
588	ndims, ndims_exp, ndims_max, &
589	lendim, lendim_exp, lendim_max )
590	if (ierr .ne. 0) then
591	write(*,9100) 'error doing ncread of PB', fnamenc
592	stop
593	end if
594
595	varname = 'PH'
596	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
597	varname, dumarray_pp(0,0,kbgn), &
598	itime, &
599	ndims, ndims_exp, ndims_max, &
600	lendim, lendim_exp, lendim_max )
601	if (ierr .ne. 0) then
602	write(*,9100) 'error doing ncread of P', fnamenc
603	stop
604	end if
605
606	do k = kbgn, nwz+add_sfc_level
607	do j = 0, nymin1
608	do i = 0, nxmin1
609	zzh(i,j,k,n) = &
610	(zzh(i,j,k,n) + dumarray_pp(i,j,k))/9.81
611	end do
612	end do
613	end do
614
615	! now use dumarray_pp to store 1/density for stress calculation below
616	! if(sfc_option .eq. sfc_option_wrf) then
617
618	! varname = 'ALT'
619	! lendim_exp(3) = nuvz-add_sfc_level
620	! lendim_max(3) = nwzmax
621	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
622	! varname, dumarray_pp(0,0,kbgn), &
623	! itime, &
624	! ndims, ndims_exp, ndims_max, &
625	! lendim, lendim_exp, lendim_max )
626	! if (ierr .ne. 0) then
627	! write(*,9100) 'error doing ncread of ALT', fnamenc
628	! stop
629	! end if
630
631	! end if
632
633	! temperature - read perturbation potential temperature,
634	! add 300. (base value), then add and convert
635	! the wrf output file contains (nuvz-add_sfc_level) levels
636	! read the data into k=kbgn,nuvz
637	varname = 'T'
638	lendim_exp(3) = nuvz-add_sfc_level
639	lendim_max(3) = nuvzmax
640	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
641	varname, tth(0,0,kbgn,n), &
642	itime, &
643	ndims, ndims_exp, ndims_max, &
644	lendim, lendim_exp, lendim_max )
645	if (ierr .ne. 0) then
646	write(*,9100) 'error doing ncread of T', fnamenc
647	stop
648	end if
649
650	do k = kbgn, nuvz
651	do j = 0, nymin1
652	do i = 0, nxmin1
653	! save potential temperature to ptth
654	ptth(i,j,k,n)=tth(i,j,k,n)+300.
655	tth(i,j,k,n) = (tth(i,j,k,n) + 300.) * &
656	(pph(i,j,k,n)/1.0e5)**0.286
657
658	end do
659	end do
660	end do
661
662	!-
663	if (turb_option .eq. turb_option_tke .or. &
664	turb_option .eq. turb_option_mytke ) then
665	print*, 'READ TKE',turb_option
666	! TKE - read Turbulent Kinetic,
667	! the wrf output file contains (nuvz-add_sfc_level) levels
668	! read the data into k=kbgn,nuvz
669	varname = 'TKE'
670	lendim_exp(3) = nuvz-add_sfc_level
671	lendim_max(3) = nuvzmax
672	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
673	varname, tkeh(0,0,kbgn,n), &
674	itime, &
675	ndims, ndims_exp, ndims_max, &
676	lendim, lendim_exp, lendim_max )
677
678	if (ierr .ne. 0) then
679	print*,'NO TKE available. Try TKE_PBL instead'
680	varname = 'TKE_PBL'
681	lendim_exp(3) = nuvz-add_sfc_level
682	lendim_max(3) = nuvzmax
683	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
684	varname, tkeh(0,0,kbgn,n), &
685	itime, &
686	ndims, ndims_exp, ndims_max, &
687	lendim, lendim_exp, lendim_max )
688	endif
689	if (ierr .ne. 0) then
690	print*,'NO TKE_PBL available. Try QKE instead'
691	varname = 'qke'
692	lendim_exp(3) = nuvz-add_sfc_level
693	lendim_max(3) = nuvzmax
694	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
695	varname, tkeh(0,0,kbgn,n), &
696	itime, &
697	ndims, ndims_exp, ndims_max, &
698	lendim, lendim_exp, lendim_max )
699	tkeh=tkeh/2. !conversion of qke
700	endif
701	if (ierr .ne. 0) then
702	varname = 'QKE'
703	lendim_exp(3) = nuvz-add_sfc_level
704	lendim_max(3) = nuvzmax
705	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
706	varname, tkeh(0,0,kbgn,n), &
707	itime, &
708	ndims, ndims_exp, ndims_max, &
709	lendim, lendim_exp, lendim_max )
710	tkeh=tkeh/2. !conversion of qke
711	endif
712
713	if (ierr .ne. 0) then
714	write(*,9100) 'error doing ncread of TKE', fnamenc
715	write(,)'change turb_option NOT to use TKE or change input file'
716	print*,'change SFC_OPTION to 0 as well'
717	stop
718	end if
719
720	endif
721
722
723	!-
724
725	! specific humidity - read mixing ratio (kg-water-vapor/kg-dry-air),
726	! then convert to (kg-water-vapor/kg-moist-air)
727	! the wrf output file contains (nuvz-add_sfc_level) levels
728	! read the data into k=kbgn,nuvz
729	varname = 'QVAPOR'
730	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
731	varname, qvh(0,0,kbgn,n), &
732	itime, &
733	ndims, ndims_exp, ndims_max, &
734	lendim, lendim_exp, lendim_max )
735	if (ierr .ne. 0) then
736	write(*,9100) 'error doing ncread of QVAPOR', fnamenc
737	stop
738	end if
739
740	do k = kbgn, nuvz
741	do j = 0, nymin1
742	do i = 0, nxmin1
743	qvh(i,j,k,n) = max( qvh(i,j,k,n), 0.0 )
744	qvh(i,j,k,n) = qvh(i,j,k,n)/(1.0 + qvh(i,j,k,n))
745	end do
746	end do
747	end do
748
749
750	! surface pressure
751	varname = 'PSFC'
752	lendim_exp(3) = 0
753	lendim_max(3) = 1
754	ndims_exp = 3
755	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
756	varname, ps(0,0,1,n), &
757	itime, &
758	ndims, ndims_exp, ndims_max, &
759	lendim, lendim_exp, lendim_max )
760	if (ierr .ne. 0) then
761	write(*,9100) 'error doing ncread of PSFC', fnamenc
762	stop
763	end if
764
765	! for the mexico city grid 3 simulation, the surface and
766	! level 1 pressures are not as consistent as one would like,
767	! with the problems occuring near the domain boundaries.
768	! so do the following
769	! -- calculate surface pressure from lowest level pressure, temp, height
770	! -- use wrf pressures (pph array) wherever possible
771	! (avoid using surface pressure and the akz/bkz, akm/bkm)
772	do j = 0, nymin1
773	do i = 0, nxmin1
774	duma = ps(i,j,1,n)
775	dumdz = 0.5*(zzh(i,j,kbgn+1,n) - zzh(i,j,kbgn,n))
776	tv = tth(i,j,kbgn,n)(1.+0.61qvh(i,j,kbgn,n))
777	ps(i,j,1,n) = pph(i,j,kbgn,n)exp( dumdzga/(r_air*tv) )
778	end do
779	end do
780
781
782	! 10 meter u velocity
783	! note: u10 is on the "T-grid" already
784	varname = 'U10'
785	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
786	varname, u10(0,0,1,n), &
787	itime, &
788	ndims, ndims_exp, ndims_max, &
789	lendim, lendim_exp, lendim_max )
790	if (ierr .ne. 0) then
791	write(*,9100) 'error doing ncread of U10', fnamenc
792	stop
793	end if
794
795
796	! 10 meter v velocity
797	! note: v10 is on the "T-grid" already
798	varname = 'V10'
799	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
800	varname, v10(0,0,1,n), &
801	itime, &
802	ndims, ndims_exp, ndims_max, &
803	lendim, lendim_exp, lendim_max )
804	if (ierr .ne. 0) then
805	write(*,9100) 'error doing ncread of V10', fnamenc
806	stop
807	end if
808
809
810	! 2 meter temperature
811	varname = 'T2'
812	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
813	varname, tt2(0,0,1,n), &
814	itime, &
815	ndims, ndims_exp, ndims_max, &
816	lendim, lendim_exp, lendim_max )
817	if (ierr .ne. 0) then
818	write(*,9100) 'error doing ncread of T2', fnamenc
819	stop
820	end if
821
822
823	! 2 meter dew point - read 2 meter water vapor mixing ratio
824	! then calculate the dew point
825	varname = 'Q2'
826	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
827	varname, td2(0,0,1,n), &
828	itime, &
829	ndims, ndims_exp, ndims_max, &
830	lendim, lendim_exp, lendim_max )
831	if (ierr .ne. 0) then
832	write(*,9100) 'error doing ncread of Q2'
833	!, fnamenc
834	do j = 0, nymin1
835	do i = 0, nxmin1
836	! 29-nov-2005 - changed qvh(i,j,1,n) to qvh(i,j,kbgn,n) here
837	td2(i,j,1,n) = qvh(i,j,kbgn,n)
838	end do
839	end do
840	end if
841
842	if (wind_option.ge.1) then
843	! print*,'mean wind from WRF is used'
844	! print*,'option ',wind_option
845	varname = 'MU'
846
847	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
848	varname, mu(0,0,1), &
849	itime, &
850	ndims, ndims_exp, ndims_max, &
851	lendim, lendim_exp, lendim_max )
852	if (ierr .ne. 0) then
853	write(*,9100) 'error doing MU', fnamenc
854	stop
855	end if
856
857	varname = 'MUB'
858
859	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
860	varname, mub(0,0,1), &
861	itime, &
862	ndims, ndims_exp, ndims_max, &
863	lendim, lendim_exp, lendim_max )
864	if (ierr .ne. 0) then
865	write(*,9100) 'error doing MUV', fnamenc
866	stop
867	end if
868	endif
869
870	! varname = 'MAPFAC_UX'
871	! varname = 'MAPFAC_UY' !try
872	! lendim_exp(1) = nx+1
873	! lendim_max(1) = nxmax
874	! lendim_exp(2) = ny
875	! lendim_max(2) = nymax
876	!
877	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
878	! varname, m_u(0,0,1), &
879	! itime, &
880	! ndims, ndims_exp, ndims_max, &
881	! lendim, lendim_exp, lendim_max )
882	! if (ierr .ne. 0) then
883	! write(*,9100) 'error doing MAP U', fnamenc
884	! stop
885	! end if
886	!
887	! varname = 'MAPFAC_VY'
888	! varname = 'MAPFAC_VX' !try
889	! lendim_exp(1) = nx
890	! lendim_max(1) = nxmax
891	! lendim_exp(2) = ny+1
892	! lendim_max(2) = nymax
893	!
894	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
895	! varname, m_v(0,0,1), &
896	! itime, &
897	! ndims, ndims_exp, ndims_max, &
898	! lendim, lendim_exp, lendim_max )
899	! if (ierr .ne. 0) then
900	! write(*,9100) 'error doing MAP V', fnamenc
901	! stop
902	! end if
903
904	! varname = 'MAPFAC_MX'
905	! lendim_exp(1) = nx
906	! lendim_max(1) = nxmax
907	! lendim_exp(2) = ny
908	! lendim_max(2) = nymax
909	!
910	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
911	! varname, m_x(0,0,1), &
912	! itime, &
913	! ndims, ndims_exp, ndims_max, &
914	! lendim, lendim_exp, lendim_max )
915	! if (ierr .ne. 0) then
916	! write(*,9100) 'error doing MAP X', fnamenc
917	! varname = 'MAPFAC_U'
918	! lendim_exp(1) = nx+1
919	! lendim_max(1) = nxmax
920	! lendim_exp(2) = ny
921	! lendim_max(2) = nymax
922	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
923	! varname, m_u(0,0,1), &
924	! itime, &
925	! ndims, ndims_exp, ndims_max, &
926	! lendim, lendim_exp, lendim_max )
927	! do j = 0, nymin1
928	! do i = 0, nxmin1
929	! m_x(i,j,1)=(m_u(i,j,1)+m_u(i+1,j,1))*0.5
930	! enddo
931	! enddo
932	! if (ierr .ne. 0) then
933	! write(*,9100) 'error doing MAP U', fnamenc
934	! print*,'NO MAP FACTOR IS GOING TO BE USED.'
935	! print*,'LARGE UNCERTAINTIES TO BE EXPECTED'
936	! do j = 0, nymin1
937	! do i = 0, nxmin1
938	! m_x(i,j,1)=1.
939	! enddo
940	! enddo
941	! end if
942	! end if
943
944	! varname = 'MAPFAC_M'
945	! lendim_exp(1) = nx
946	! lendim_max(1) = nxmax
947	! lendim_exp(2) = ny
948	! lendim_max(2) = nymax
949	!
950	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
951	! varname, m_z(0,0,1), &
952	! itime, &
953	! ndims, ndims_exp, ndims_max, &
954	! lendim, lendim_exp, lendim_max )
955	! if (ierr .ne. 0) then
956	! write(*,9100) 'error doing MAP W', fnamenc
957	! stop
958	! end if
959
960	! varname = 'MAPFAC_MY'
961	! lendim_exp(1) = nx
962	! lendim_max(1) = nxmax
963	! lendim_exp(2) = ny
964	! lendim_max(2) = nymax
965	!
966	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
967	! varname, m_y(0,0,1), &
968	! itime, &
969	! ndims, ndims_exp, ndims_max, &
970	! lendim, lendim_exp, lendim_max )
971	! if (ierr .ne. 0) then
972	! write(*,9100) 'error doing MAP Y', fnamenc
973	! varname = 'MAPFAC_V'
974	! lendim_exp(1) = nx
975	! lendim_max(1) = nxmax
976	! lendim_exp(2) = ny+1
977	! lendim_max(2) = nymax
978	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
979	! varname, m_v(0,0,1), &
980	! itime, &
981	! ndims, ndims_exp, ndims_max, &
982	! lendim, lendim_exp, lendim_max )
983	! do j = 0, nymin1
984	! do i = 0, nxmin1
985	! m_y(i,j,1)=(m_v(i,j,1)+m_v(i,j+1,1))*0.5
986	! enddo
987	! enddo
988	! if (ierr .ne. 0) then
989	! write(*,9100) 'ERROR doing MAP V', fnamenc
990	! print*,'NO MAP FACTOR IS GOING TO BE USED.'
991	! print*,'LARGE UNCERTAINTIES TO BE EXPECTED'
992	! do j = 0, nymin1
993	! do i = 0, nxmin1
994	! m_y(i,j,1)=1.
995	! enddo
996	! enddo
997	! end if
998	! end if
999	lendim_exp(1) = nx
1000	lendim_max(1) = nxmax
1001	lendim_exp(2) = ny
1002	lendim_max(2) = nymax
1003
1004
1005
1006	! calculate water vapor pressure in mb, from sfc pressure
1007	! and 2 m mixing ratio
1008	iduma = 0
1009	do j = 0, nymin1
1010	do i = 0, nxmin1
1011	! 29-nov-2005 - added this to catch occasional tt2n=0.0 values
1012	duma = max( 100.0, tth(i,j,kbgn,n)-50.0 )
1013	if (tt2(i,j,1,n) .le. duma) then
1014	iduma = iduma + 1
1015	if (iduma .eq. 1) then
1016	write(,) 'readwind - bad tt2 at'
1017	write(,) 'i, j, tt2 =', i, j, tt2(i,j,1,n)
1018	end if
1019	! stop
1020	tt2(i,j,1,n) = tth(i,j,kbgn,n)
1021	td2(i,j,1,n) = qvh(i,j,kbgn,n)
1022	end if
1023	duma = td2(i,j,1,n)/0.622
1024	ewater_mb = 0.01( 0.99976ps(i,j,1,n)*duma/(1.0+duma) )
1025	esatwater_mb = 0.01*ew(tt2(i,j,1,n))
1026	ewater_mb = max( 1.0e-10, min( esatwater_mb, ewater_mb ) )
1027	! then use the following, which is from an old 1970's report
1028	! (reference not available, but the formula works)
1029	! tdew(in C) = (4318.76/(19.5166 - ln(ewater(in mb)))) - 243.893
1030	td2(i,j,1,n) = 273.16 + &
1031	(4318.76/(19.5166 - log(ewater_mb))) - 243.893
1032	end do
1033	end do
1034	if (iduma .gt. 0) write(,) &
1035	'readwind - bad tt2 count =', iduma
1036
1037
1038	! sea level pressure - calculate it from surface pressure and
1039	! ground elevation using standard atmosphere relations
1040	do j = 0, nymin1
1041	do i = 0, nxmin1
1042	msl(i,j,1,n) = ps(i,j,1,n)/ &
1043	((1.0 - 6.5e-3oro(i,j)/288.0)*5.2553)
1044	end do
1045	end do
1046
1047
1048	! large scale precipitation
1049	! convective precipitation
1050	! the wrf output files contain these as "accumulated totals"
1051	! I need to find out if these are accumulated over the output
1052	! file frequency, or over the total run.
1053	! For now, set to zero
1054	! total cloud cover
1055	! Doesn't appear to be any 2-d cloud cover field in the
1056	! wrf output.
1057	! For now, set to zero
1058	do j = 0, nymin1
1059	do i = 0, nxmin1
1060	lsprec(i,j,1,n) = 0.0
1061	convprec(i,j,1,n) = 0.0
1062	tcc(i,j,1,n) = 0.0
1063	end do
1064	end do
1065
1066	!
1067	! Large-scale precipitation (accumulated rain, mm)
1068	! will convert to mm/h in interpolat_rain.f
1069	varname = 'RAINNC'
1070	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1071	varname, lsprec(0,0,1,n), &
1072	itime, &
1073	dims, ndims_exp, ndims_max, &
1074	lendim, lendim_exp, lendim_max )
1075	if (ierr .ne. 0) then
1076	write(*,9100) 'error doing ncread of RAINNC,set to zero', fnamenc
1077	do j = 0, nymin1
1078	do i = 0, nxmin1
1079	lsprec(i,j,1,n) = 0.0
1080	end do
1081	end do
1082	end if
1083
1084	!
1085	! Convective cumulus precipitation (accumulated rain, mm)
1086
1087	varname = 'RAINC'
1088	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1089	varname, convprec(0,0,1,n), &
1090	itime, &
1091	ndims, ndims_exp, ndims_max, &
1092	lendim, lendim_exp, lendim_max )
1093	if (ierr .ne. 0) then
1094	write(*,9100) 'error doing ncread of RAINC, set to zero', fnamenc
1095	do j = 0, nymin1
1096	do i = 0, nxmin1
1097	convprec(i,j,1,n) = 0.0
1098	end do
1099	end do
1100	end if
1101
1102	!
1103	! Clound fraction (cloud cover)
1104	! varname = 'CLDFRA'
1105	! call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1106	! varname, tcc(0,0,1,n), &
1107	! itime, &
1108	! ndims, ndims_exp, ndims_max, &
1109	! lendim, lendim_exp, lendim_max )
1110	! if (ierr .ne. 0) then
1111	!! write(*,9100) 'error doing ncread of CLDFRA, set to zero'
1112	!!, fnamenc
1113	! do j = 0, nymin1
1114	! do i = 0, nxmin1
1115	! tcc(i,j,1,n) = 0.0
1116	! end do
1117	! end do
1118	! end if
1119	!!C write(,)'read CLDFRA 0-sucess ',ierr
1120
1121	!
1122
1123	! snow depth
1124	varname = 'SNOWH'
1125	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1126	varname, sd(0,0,1,n), &
1127	itime, &
1128	ndims, ndims_exp, ndims_max, &
1129	lendim, lendim_exp, lendim_max )
1130	if (ierr .ne. 0) then
1131	! write(*,9100) 'error doing ncread of SNOWH', fnamenc
1132	do j = 0, nymin1
1133	do i = 0, nxmin1
1134	sd(i,j,1,n) = 0.0
1135	end do
1136	end do
1137	end if
1138
1139
1140	! surface sensible heat flux (positive <--> upwards)
1141	varname = 'HFX'
1142	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1143	varname, sshf(0,0,1,n), &
1144	itime, &
1145	ndims, ndims_exp, ndims_max, &
1146	lendim, lendim_exp, lendim_max )
1147	do j = 0, nymin1
1148	do i = 0, nxmin1
1149	sshf(i,j,1,n) = -sshf(i,j,1,n)
1150	end do
1151	end do
1152
1153	if (ierr .ne. 0) then
1154	write(*,9100) 'error doing ncread of HFX', fnamenc
1155	do j = 0, nymin1
1156	do i = 0, nxmin1
1157	sshf(i,j,1,n) = 0.0
1158	end do
1159	end do
1160	hflswitch=.false. ! Heat flux is not available
1161	else
1162	hflswitch=.true. ! Heat flux is available
1163	! limit to values to bounds originally used by flexpart?
1164	! do 1502 j=0,nymin1
1165	! do 1502 i=0,nxmin1
1166	! if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
1167	! if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
1168	!1502 continue
1169	end if
1170
1171	! ustar
1172	varname = 'UST'
1173	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1174	varname, ustar(0,0,1,n), &
1175	itime, &
1176	ndims, ndims_exp, ndims_max, &
1177	lendim, lendim_exp, lendim_max )
1178	if (ierr .ne. 0) then
1179	write(*,9100) 'error doing ncread of UST', fnamenc
1180	do j = 0, nymin1
1181	do i = 0, nxmin1
1182	ustar(i,j,1,n) = 0.0
1183	end do
1184	end do
1185	strswitch=.false. ! ustar is not available
1186	else
1187	strswitch=.true. ! ustar is available
1188	do j=0,nymin1
1189	do i=0,nxmin1
1190	! surfstr(i,j,1,n)=ustar(i,j,1,n)/dumarray_pp(i,j,kbgn)
1191	enddo
1192	enddo
1193	end if
1194
1195	! pblh
1196	if(sfc_option .eq. sfc_option_wrf) then
1197	varname = 'PBLH'
1198	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1199	varname, hmix(0,0,1,n), &
1200	itime, &
1201	ndims, ndims_exp, ndims_max, &
1202	lendim, lendim_exp, lendim_max )
1203	if (ierr .ne. 0) then
1204	write(*,9100) 'error doing ncread of PBLH', fnamenc
1205	stop
1206	endif
1207
1208	endif
1209
1210	! surface solar radiation flux (positive <--> downwards)
1211	varname = 'SWDOWN'
1212	call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
1213	varname, ssr(0,0,1,n), &
1214	itime, &
1215	ndims, ndims_exp, ndims_max, &
1216	lendim, lendim_exp, lendim_max )
1217	if (ierr .ne. 0) then
1218	write(*,9100) 'error doing ncread of SWDOWN', fnamenc
1219	do j = 0, nymin1
1220	do i = 0, nxmin1
1221	ssr(i,j,1,n) = 0.0
1222	end do
1223	end do
1224	else
1225	do j = 0, nymin1
1226	do i = 0, nxmin1
1227	ssr(i,j,1,n) = max( ssr(i,j,1,n), 0.0 )
1228	end do
1229	end do
1230	end if
1231
1232
1233	! ew & ns surface stress
1234	! Doesn't appear to be any 2-d cloud cover field in the
1235	! wrf output.
1236	! For now, set to zero
1237	do j = 0, nymin1
1238	do i = 0, nxmin1
1239	ewss(i,j) = 0.0
1240	nsss(i,j) = 0.0
1241	end do
1242	end do
1243	! strswitch=.false. ! Surface stress is not available
1244
1245
1246	! orography
1247	! standard deviation of orography
1248	! land sea mask
1249	! these should be fixed during a simulation
1250	! so there is no reason to do them again ??
1251
1252
1253	! * done with reading the wrf output file *
1254
1255	! print*,'uu out1',uuh(0,259,1:10)
1256	! print*,'mu out1',mu(0,259,1),mub(0,259,1)
1257	! print*,'m_xn out1',m_x(0,259,1),m_y(0,259,1)
1258
1259
1260	! interpolate uuh from the "U-grid" to the "T-grid"
1261	! interpolate vvh from the "V-grid" to the "T-grid"
1262	! new: convert mass weighted wind to wind.
1263	! print*,'wind_option',wind_option
1264	if (wind_option.le.0) then
1265	! if (wind_option.eq.-1) then
1266	! call calc_uvmet(uuh,vvh,urot,vrot,1)
1267	! endif
1268	do k = kbgn, nuvz
1269	do j = 0, nymin1
1270	do i = 0, nxmin1
1271	! needs to rotate u and v to work. needs read alpha or something like that.
1272	! if in mercator, no need.
1273	! divh(i,j,k)=(uuh(i+1,j,k)-uuh(i,j,k))/dx &
1274	! +(vvh(i,j+1,k)-vvh(i,j,k))/dy
1275	if (wind_option.lt.0) then
1276	divh(i,j,k)=(uuh(i+1,j,k)-uuh(i,j,k))/dx*m_x(i,j,1) &
1277	+(vvh(i,j+1,k)-vvh(i,j,k))/dy*m_y(i,j,1)
1278	endif
1279	uuh(i,j,k) = 0.5*(uuh(i,j,k) + uuh(i+1,j,k))
1280	vvh(i,j,k) = 0.5*(vvh(i,j,k) + vvh(i,j+1,k))
1281	end do
1282	end do
1283	end do
1284	elseif (wind_option.eq.1) then
1285	do k = kbgn, nuvz
1286	do j = 0, nymin1
1287	do i = 0, nxmin1
1288	! uuh(i,j,k) = 0.5(uuh(i,j,k)m_u(i,j,1) + uuh(i+1,j,k)*m_u(i+1,j,1))
1289	! vvh(i,j,k) = 0.5(vvh(i,j,k)m_v(i,j,1) + vvh(i,j+1,k)*m_v(i,j+1,1))
1290	uuh(i,j,k) = 0.5*(uuh(i,j,k) + uuh(i+1,j,k))
1291	vvh(i,j,k) = 0.5*(vvh(i,j,k) + vvh(i,j+1,k))
1292	mu2=mu(i,j,1)+mub(i,j,1)
1293	! uuh(i,j,k) = uuh(i,j,k)/mu2!0.5(m_u(i,j,1)+m_u(i+1,j,1))
1294	! vvh(i,j,k) = vvh(i,j,k)/mu2!0.5(m_v(i,j,1)+m_v(i,j+1,1))
1295	! m_y=0.5*(m_v(i,j,1)+m_v(i,j+1,1))
1296	uuh(i,j,k) = uuh(i,j,k)/mu2 !*m_y(i,j,1) !without m=true wind
1297	vvh(i,j,k) = vvh(i,j,k)/mu2 !*m_x(i,j,1)
1298	wwh(i,j,k) = wwh(i,j,k)/mu2 !*m_y(i,j,1)
1299	end do
1300	end do
1301	end do
1302	elseif (wind_option.eq.2) then
1303	do k = kbgn, nuvz
1304	do j = 0, nymin1
1305	do i = 0, nxmin1
1306	uuh(i,j,k) = 0.5*(uuh(i,j,k) + uuh(i+1,j,k))
1307	vvh(i,j,k) = 0.5*(vvh(i,j,k) + vvh(i,j+1,k))
1308	mu2=mu(i,j,1)+mub(i,j,1)
1309	wwh(i,j,k) = wwh(i,j,k)/mu2 !*m_y(i,j,1)
1310	end do
1311	end do
1312	end do
1313
1314	endif
1315
1316
1317	! for ecmwf flexpart, if nwz = nlev_ec+1, then wwh is set
1318	! to zero at the top level
1319	! for wrf, nlev_ec==n_bottom_top and so nwz = nlev_ec+1.
1320	! however, it doesn't seem appropriate to zero wwh at
1321	! the model top which might be ~100 hPa.
1322	! so deactivate this for now
1323	! if(levdiff2.eq.0) then
1324	! iwmax=nlev_ec+1
1325	! do 60 i=0,nxmin1
1326	! do 60 j=0,nymin1
1327	!60 wwh(i,j,nlev_ec+1)=0.
1328	! endif
1329
1330	! For global fields, assign the leftmost data column also to the rightmost
1331	! data column; if required, shift whole grid by nxshift grid points
1332	!
1333	! FLEXPART_WRF - all "global" stuff is turned off
1334	!*************************************************************************
1335
1336	! if (xglobal) then
1337	! call shift_field_0(ewss,nxfield,ny)
1338	! call shift_field_0(nsss,nxfield,ny)
1339	! call shift_field_0(oro,nxfield,ny)
1340	! call shift_field_0(excessoro,nxfield,ny)
1341	! call shift_field_0(lsm,nxfield,ny)
1342	! call shift_field(ps,nxfield,ny,1,1,2,n)
1343	! call shift_field(sd,nxfield,ny,1,1,2,n)
1344	! call shift_field(msl,nxfield,ny,1,1,2,n)
1345	! call shift_field(tcc,nxfield,ny,1,1,2,n)
1346	! call shift_field(u10,nxfield,ny,1,1,2,n)
1347	! call shift_field(v10,nxfield,ny,1,1,2,n)
1348	! call shift_field(tt2,nxfield,ny,1,1,2,n)
1349	! call shift_field(td2,nxfield,ny,1,1,2,n)
1350	! call shift_field(lsprec,nxfield,ny,1,1,2,n)
1351	! call shift_field(convprec,nxfield,ny,1,1,2,n)
1352	! call shift_field(sshf,nxfield,ny,1,1,2,n)
1353	! call shift_field(ssr,nxfield,ny,1,1,2,n)
1354	! call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
1355	! call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
1356	! call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
1357	! call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
1358	! call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
1359	! endif
1360
1361	! CALCULATE SURFSTR
1362	if(sfc_option .eq. sfc_option_diagnosed) then
1363	do i=0,nxmin1
1364	do j=0,nymin1
1365	surfstr(i,j,1,n)=sqrt(ewss(i,j)2+nsss(i,j)2)
1366	enddo
1367	enddo
1368	strswitch=.false.
1369	endif
1370
1371	if ((.not.hflswitch).or.(.not.strswitch)) then
1372	write(,) 'WARNING: No (or incomplete) flux data ' // &
1373	'contained in WRF output file ', &
1374	wfname(indj)
1375
1376	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
1377	! As ECMWF has increased the model resolution, such that now the first model
1378	! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
1379	! (3rd model level in FLEXPART) for the profile method
1380	!
1381	! FLEXPART_WRF - use k=(2+add_sfc_level) here instead of k=3
1382	!***************************************************************************
1383
1384	k = 2 + add_sfc_level
1385	do i=0,nxmin1
1386	do j=0,nymin1
1387	! plev1=akz(3)+bkz(3)*ps(i,j,1,n)
1388	plev1=pph(i,j,k,n)
1389	pmean=0.5*(ps(i,j,1,n)+plev1)
1390	tv=tth(i,j,k,n)(1.+0.61qvh(i,j,k,n))
1391	fu=-r_air*tv/ga/pmean
1392	hlev1=fu*(plev1-ps(i,j,1,n)) ! HEIGTH OF FIRST MODEL LAYER
1393	ff10m= sqrt(u10(i,j,1,n)2+v10(i,j,1,n)2)
1394	fflev1=sqrt(uuh(i,j,k)2+vvh(i,j,k)2)
1395	call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
1396	tt2(i,j,1,n),tth(i,j,k,n),ff10m,fflev1, &
1397	surfstr(i,j,1,n),sshf(i,j,1,n))
1398	if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
1399	if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
1400	enddo
1401	enddo
1402	endif
1403
1404
1405	! Assign 10 m wind to model level at eta=1.0 to have one additional model
1406	! level at the ground
1407	! Specific humidity is taken the same as at one level above
1408	! Temperature is taken as 2 m temperature
1409	!
1410	! Note that the uuh, vvh, tth, & qvh data have already been shifted
1411	! upwards by one level, when they were read in.
1412	!**************************************************************************
1413
1414	if (add_sfc_level .eq. 1) then
1415	do j = 0, nymin1
1416	do i = 0, nxmin1
1417	uuh(i,j,1) = u10(i,j,1,n)
1418	vvh(i,j,1) = v10(i,j,1,n)
1419	tth(i,j,1,n) = tt2(i,j,1,n)
1420	ptth(i,j,1,n) = ptth(i,j,2,n)
1421	qvh(i,j,1,n) = qvh(i,j,2,n)
1422	tkeh(i,j,1,n)=tkeh(i,j,2,n)
1423	! pressure at 2 m AGL
1424	pph(i,j,1,n) = 0.99976*ps(i,j,1,n)
1425	! height (MSL) at ground level (shift it down)
1426	zzh(i,j,1,n) = zzh(i,j,2,n)
1427	! height (MSL) at top of the added level
1428	zzh(i,j,2,n) = zzh(i,j,1,n) + 4.0
1429	if (hmix(i,j,1,n).lt.hmixmin) hmix(i,j,1,n)=hmixmin
1430
1431	enddo
1432	enddo
1433	end if
1434
1435
1436	do i=0,nxmax-1
1437	do j=0,nymax-1
1438	do k=1,nuvzmax
1439	u_wrf(i,j,k,n)=uuh(i,j,k)
1440	v_wrf(i,j,k,n)=vvh(i,j,k)
1441	enddo
1442	enddo
1443	enddo
1444
1445	do i=0,nxmax-1
1446	do j=0,nymax-1
1447	do k=1,nwzmax
1448	w_wrf(i,j,k,n)=wwh(i,j,k)
1449	enddo
1450	enddo
1451	enddo
1452
1453
1454
1455
1456
1457	return
1458	end subroutine readwind
1459

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

Download in other formats:

Original Format