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.61*qvhn(i,j,kbgn,n,l)) |
---|
731 | psn(i,j,1,n,l) = pphn(i,j,kbgn,n,l)*exp( dumdz*ga/(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.99976*psn(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-3*oron(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.61*qvhn(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 |
---|