Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

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

Download in other formats:

Original Format