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.61*qvh(i,j,kbgn,n)) |
---|
777 | ps(i,j,1,n) = pph(i,j,kbgn,n)*exp( dumdz*ga/(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.99976*ps(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-3*oro(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.61*qvh(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 | |
---|