1 | subroutine gridcheck_nests |
---|
2 | !******************************************************************************* |
---|
3 | ! * |
---|
4 | ! This routine checks the grid specification for the nested model domains. * |
---|
5 | ! It is similar to subroutine gridcheck, which checks the mother domain. * |
---|
6 | ! * |
---|
7 | ! Note: This is the FLEXPART_WRF version of subroutine gridcheck. * |
---|
8 | ! The computational grid is the WRF x-y grid rather than lat-lon. * |
---|
9 | ! There are many differences from the FLEXPART version. * |
---|
10 | ! * |
---|
11 | ! Authors: A. Stohl, G. Wotawa * |
---|
12 | ! 8 February 1999 * |
---|
13 | ! * |
---|
14 | ! Nov 2005 - R. Easter - MAJOR revisions for FLEXPART_WRF * |
---|
15 | ! Dec 2005, R. Easter - changed names of "*lon0*" & "*lat0*" variables * |
---|
16 | ! * |
---|
17 | !******************************************************************************* |
---|
18 | |
---|
19 | ! use grib_api |
---|
20 | use par_mod |
---|
21 | use com_mod |
---|
22 | |
---|
23 | |
---|
24 | integer,parameter :: ndims_max=4 |
---|
25 | integer :: i, ierr, ifn, itime, ix |
---|
26 | integer :: idiagaa, idiagaa_1, idiagaa_2, idiagbb |
---|
27 | integer :: iduma, idumb |
---|
28 | integer :: jy |
---|
29 | integer :: k |
---|
30 | integer :: l, lp |
---|
31 | integer :: lendim(ndims_max), lendim_exp(ndims_max), & |
---|
32 | lendim_max(ndims_max) |
---|
33 | integer :: m |
---|
34 | integer :: map_proj_id_dum |
---|
35 | integer :: ndims, ndims_exp, & |
---|
36 | ext_scalar,pbl_physics,mp_physics_dum |
---|
37 | integer :: n_west_east, n_south_north, n_bottom_top |
---|
38 | integer :: nuvzn, nwzn |
---|
39 | |
---|
40 | real :: dx_met, dy_met |
---|
41 | real :: duma, dumb, dumc, dumx, dumy |
---|
42 | real :: dump1, dump2, dumdz |
---|
43 | real :: eta_w_wrf_nest(nwzmax), eta_u_wrf_nest(nwzmax) |
---|
44 | real :: map_stdlon_dum, map_truelat1_dum, map_truelat2_dum |
---|
45 | real :: pint, p_top_wrf_nest |
---|
46 | real :: xaux1, xaux2, yaux1, yaux2 |
---|
47 | |
---|
48 | character(len=160) :: fnamenc, varname |
---|
49 | |
---|
50 | real(kind=4) :: m_un(0:nxmaxn-1,0:nymaxn-1,1,maxnests) |
---|
51 | real(kind=4) :: m_vn(0:nxmaxn-1,0:nymaxn-1,1,maxnests) |
---|
52 | |
---|
53 | ! Loop about all nesting levels |
---|
54 | !****************************** |
---|
55 | ! idiagaa_1 = 1 |
---|
56 | idiagaa_1 = 0 |
---|
57 | idiagaa_2 = 0 |
---|
58 | ! idiagbb = 10 |
---|
59 | idiagbb = 0 |
---|
60 | |
---|
61 | do l=1,numbnests |
---|
62 | |
---|
63 | write(*,'(//a,i3)') 'gridcheck_nests output for grid #', l |
---|
64 | |
---|
65 | ! |
---|
66 | ! get grid info from the wrf netcdf file |
---|
67 | ! |
---|
68 | if(ideltas.gt.0) then |
---|
69 | ifn=1 |
---|
70 | else |
---|
71 | ifn=numbwf |
---|
72 | endif |
---|
73 | m = numpath+2*(l-1)+1 |
---|
74 | fnamenc = path(m)(1:length(m)) // wfnamen(l,ifn) |
---|
75 | |
---|
76 | idiagaa = idiagaa_1 |
---|
77 | |
---|
78 | call read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, & |
---|
79 | n_west_east, n_south_north, n_bottom_top, & |
---|
80 | dx_met, dy_met, & |
---|
81 | m_grid_id(l), m_parent_grid_id(l), m_parent_grid_ratio(l), & |
---|
82 | i_parent_start(l), j_parent_start(l), & |
---|
83 | map_proj_id_dum, map_stdlon_dum, & |
---|
84 | map_truelat1_dum, map_truelat2_dum, & |
---|
85 | ext_scalar,pbl_physics,mp_physics_dum ) |
---|
86 | if (ierr .ne. 0) goto 999 |
---|
87 | |
---|
88 | |
---|
89 | mp_physicsn(l)=mp_physics_dum |
---|
90 | |
---|
91 | ! subtract 1 because i & j indexing in flexpart always starts at 0 |
---|
92 | i_parent_start(l) = i_parent_start(l)-1 |
---|
93 | j_parent_start(l) = j_parent_start(l)-1 |
---|
94 | |
---|
95 | ! |
---|
96 | ! set grid dimension and size variables |
---|
97 | ! |
---|
98 | nxn(l) = n_west_east |
---|
99 | nyn(l) = n_south_north |
---|
100 | |
---|
101 | nuvzn = n_bottom_top |
---|
102 | nwzn = n_bottom_top + 1 |
---|
103 | |
---|
104 | ! for FLEXPART_WRF, x & y coords are in meters |
---|
105 | dxn(l) = dx_met |
---|
106 | dyn(l) = dy_met |
---|
107 | |
---|
108 | |
---|
109 | ! |
---|
110 | ! check that grid dimensions are not too big |
---|
111 | ! |
---|
112 | ! flexpart_wrf 07-nov-2005 - require (nxn+1 .le. nxmaxn) and (nyn+1 .le. nymaxn) |
---|
113 | ! because u,v in met. files are on staggered grid |
---|
114 | if (nxn(l)+1 .gt. nxmaxn) then |
---|
115 | write(*,*) 'FLEXPART gridcheck_nests error: ' // & |
---|
116 | 'Too many grid points in x direction.' |
---|
117 | write(*,*) 'Change parameter settings in file par_mod.' |
---|
118 | write(*,*) 'l, nxn(l)+1, nxmaxn =', l, nxn(l), nxmaxn |
---|
119 | stop |
---|
120 | endif |
---|
121 | |
---|
122 | if (nyn(l)+1 .gt. nymaxn) then |
---|
123 | write(*,*) 'FLEXPART gridcheck_nests error: ' // & |
---|
124 | 'Too many grid points in y direction.' |
---|
125 | write(*,*) 'Change parameter settings in file par_mod.' |
---|
126 | write(*,*) 'l, nyn(l)+1, nymaxn =', l, nyn(l), nymaxn |
---|
127 | stop |
---|
128 | endif |
---|
129 | |
---|
130 | nuvzn = nuvzn+add_sfc_level |
---|
131 | if (nuvzn .ne. nuvz) then |
---|
132 | write(*,*) 'FLEXPART gridcheck_nests error: ' // & |
---|
133 | 'nuvzn and nuvz differ' |
---|
134 | write(*,*) 'l, nuvzn, nuvz =', l, nuvzn, nuvz |
---|
135 | stop |
---|
136 | endif |
---|
137 | |
---|
138 | if (nwzn .ne. nwz) then |
---|
139 | write(*,*) 'FLEXPART gridcheck_nests error: ' // & |
---|
140 | 'nwzn and nwz differ' |
---|
141 | write(*,*) 'l, nwzn, nwz =', l, nwzn, nwz |
---|
142 | stop |
---|
143 | endif |
---|
144 | |
---|
145 | ! check that map projection info matches parent |
---|
146 | duma = 3.0e-7*max( abs(map_stdlon), 1.0e-30 ) |
---|
147 | dumb = 3.0e-7*max( abs(map_truelat1), 1.0e-30 ) |
---|
148 | dumc = 3.0e-7*max( abs(map_truelat2), 1.0e-30 ) |
---|
149 | iduma = 0 |
---|
150 | if (map_proj_id .ne. map_proj_id_dum) iduma = 1 |
---|
151 | if (abs(map_stdlon-map_stdlon_dum) .gt. duma) iduma = 2 |
---|
152 | if (abs(map_truelat1-map_truelat1_dum) .gt. dumb) iduma = 3 |
---|
153 | if (abs(map_truelat2-map_truelat2_dum) .gt. dumc) iduma = 4 |
---|
154 | if (iduma .ne. 0) then |
---|
155 | write(*,*) 'FLEXPART gridcheck_nests error: ' // & |
---|
156 | 'map projection parameters differ' |
---|
157 | write(*,*) 'l, map param #=', l, iduma |
---|
158 | stop |
---|
159 | end if |
---|
160 | |
---|
161 | varname = 'MAPFAC_MX' |
---|
162 | lendim_exp(1) = nxn(l) |
---|
163 | lendim_max(1) = nxmaxn |
---|
164 | lendim_exp(2) = nyn(l) |
---|
165 | lendim_max(2) = nymaxn |
---|
166 | ndims_exp = 3 |
---|
167 | itime=1 |
---|
168 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
169 | varname, m_xn(0,0,1,l), & |
---|
170 | itime, & |
---|
171 | ndims, ndims_exp, ndims_max, & |
---|
172 | lendim, lendim_exp, lendim_max ) |
---|
173 | if (ierr .ne. 0) then |
---|
174 | varname = 'MAPFAC_M' |
---|
175 | lendim_exp(1) = nxn(l) |
---|
176 | lendim_max(1) = nxmaxn |
---|
177 | lendim_exp(2) = nyn(l) |
---|
178 | lendim_max(2) = nymaxn |
---|
179 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
180 | varname, m_xn(0,0,1,l), & |
---|
181 | itime, & |
---|
182 | ndims, ndims_exp, ndims_max, & |
---|
183 | lendim, lendim_exp, lendim_max ) |
---|
184 | endif |
---|
185 | if (ierr .ne. 0) then |
---|
186 | print*,'error doing MAP X' |
---|
187 | varname = 'MAPFAC_UX' |
---|
188 | lendim_exp(1) = nxn(l)+1 |
---|
189 | lendim_max(1) = nxmaxn |
---|
190 | lendim_exp(2) = nyn(l) |
---|
191 | lendim_max(2) = nymaxn |
---|
192 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
193 | varname, m_un(0,0,1,l), & |
---|
194 | itime, & |
---|
195 | ndims, ndims_exp, ndims_max, & |
---|
196 | lendim, lendim_exp, lendim_max ) |
---|
197 | do j = 0, nyn(l)-1 |
---|
198 | do i = 0, nxn(l)-1 |
---|
199 | m_xn(i,j,1,l)=(m_un(i,j,1,l)+m_un(i+1,j,1,l))*0.5 |
---|
200 | enddo |
---|
201 | enddo |
---|
202 | if (ierr .ne. 0) then |
---|
203 | print*,'error doing MAP U' |
---|
204 | print*,'NO MAP FACTOR IS GOING TO BE USED.' |
---|
205 | print*,'LARGE UNCERTAINTIES TO BE EXPECTED' |
---|
206 | do j = 0, nyn(l)-1 |
---|
207 | do i = 0, nxn(l)-1 |
---|
208 | m_xn(i,j,1,l)=1. |
---|
209 | enddo |
---|
210 | enddo |
---|
211 | end if |
---|
212 | end if |
---|
213 | |
---|
214 | varname = 'MAPFAC_MY' |
---|
215 | lendim_exp(1) = nxn(l) |
---|
216 | lendim_max(1) = nxmaxn |
---|
217 | lendim_exp(2) = nyn(l) |
---|
218 | lendim_max(2) = nymaxn |
---|
219 | |
---|
220 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
221 | varname, m_yn(0,0,1,l), & |
---|
222 | itime, & |
---|
223 | ndims, ndims_exp, ndims_max, & |
---|
224 | lendim, lendim_exp, lendim_max ) |
---|
225 | if (ierr .ne. 0) then |
---|
226 | varname = 'MAPFAC_M' |
---|
227 | lendim_exp(1) = nxn(l) |
---|
228 | lendim_max(1) = nxmaxn |
---|
229 | lendim_exp(2) = nyn(l) |
---|
230 | lendim_max(2) = nymaxn |
---|
231 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
232 | varname, m_yn(0,0,1,l), & |
---|
233 | itime, & |
---|
234 | ndims, ndims_exp, ndims_max, & |
---|
235 | lendim, lendim_exp, lendim_max ) |
---|
236 | endif |
---|
237 | if (ierr .ne. 0) then |
---|
238 | print*,'error doing MAP Y' |
---|
239 | varname = 'MAPFAC_VY' |
---|
240 | lendim_exp(1) = nxn(l) |
---|
241 | lendim_max(1) = nxmaxn |
---|
242 | lendim_exp(2) = nyn(l)+1 |
---|
243 | lendim_max(2) = nymaxn |
---|
244 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
245 | varname, m_vn(0,0,1,l), & |
---|
246 | itime, & |
---|
247 | ndims, ndims_exp, ndims_max, & |
---|
248 | lendim, lendim_exp, lendim_max ) |
---|
249 | do j = 0, nyn(l)-1 |
---|
250 | do i = 0, nxn(l)-1 |
---|
251 | m_yn(i,j,1,l)=(m_vn(i,j,1,l)+m_vn(i,j+1,1,l))*0.5 |
---|
252 | enddo |
---|
253 | enddo |
---|
254 | if (ierr .ne. 0) then |
---|
255 | print*,'ERROR doing MAP V' |
---|
256 | print*,'NO MAP FACTOR IS GOING TO BE USED.' |
---|
257 | print*,'LARGE UNCERTAINTIES TO BE EXPECTED' |
---|
258 | do j = 0, nyn(l)-1 |
---|
259 | do i = 0, nxn(l)-1 |
---|
260 | m_yn(i,j,1,l)=1. |
---|
261 | enddo |
---|
262 | enddo |
---|
263 | end if |
---|
264 | end if |
---|
265 | lendim_exp(1) = nxn(l) |
---|
266 | lendim_max(1) = nxmaxn |
---|
267 | lendim_exp(2) = nyn(l) |
---|
268 | lendim_max(2) = nymaxn |
---|
269 | |
---|
270 | ! |
---|
271 | ! read latitude and longitude |
---|
272 | ! read oro, lsm, and excessoro |
---|
273 | ! |
---|
274 | |
---|
275 | idiagaa = idiagaa_2 |
---|
276 | |
---|
277 | varname = 'XLAT' |
---|
278 | do i = 1, ndims_max |
---|
279 | lendim_exp(i) = 0 |
---|
280 | lendim_max(i) = 1 |
---|
281 | end do |
---|
282 | itime = 1 |
---|
283 | lendim_exp(1) = nxn(l) |
---|
284 | lendim_max(1) = nxmaxn |
---|
285 | lendim_exp(2) = nyn(l) |
---|
286 | lendim_max(2) = nymaxn |
---|
287 | ndims_exp = 3 |
---|
288 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
289 | varname, ylat2dn(0,0,l), & |
---|
290 | itime, & |
---|
291 | ndims, ndims_exp, ndims_max, & |
---|
292 | lendim, lendim_exp, lendim_max ) |
---|
293 | if (ierr .ne. 0) then |
---|
294 | write(*,*) |
---|
295 | write(*,*) '*** checkgrid -- error doing ncread of XLAT' |
---|
296 | stop |
---|
297 | end if |
---|
298 | |
---|
299 | varname = 'XLONG' |
---|
300 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
301 | varname, xlon2dn(0,0,l), & |
---|
302 | itime, & |
---|
303 | ndims, ndims_exp, ndims_max, & |
---|
304 | lendim, lendim_exp, lendim_max ) |
---|
305 | if (ierr .ne. 0) then |
---|
306 | write(*,*) |
---|
307 | write(*,*) '*** checkgrid -- error doing ncread of XLONG' |
---|
308 | stop |
---|
309 | end if |
---|
310 | |
---|
311 | varname = 'HGT' |
---|
312 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
313 | varname, oron(0,0,l), & |
---|
314 | itime, & |
---|
315 | ndims, ndims_exp, ndims_max, & |
---|
316 | lendim, lendim_exp, lendim_max ) |
---|
317 | if (ierr .ne. 0) then |
---|
318 | write(*,*) |
---|
319 | write(*,*) '*** checkgrid -- error doing ncread of HGT' |
---|
320 | stop |
---|
321 | end if |
---|
322 | |
---|
323 | ! lsm = landsea mask == land fraction (or non-ocean fraction) |
---|
324 | ! for now, set lsm=1.0 which means land |
---|
325 | do jy=0,nyn(l)-1 |
---|
326 | do ix=0,nxn(l)-1 |
---|
327 | lsmn(ix,jy,l)=1.0 |
---|
328 | end do |
---|
329 | end do |
---|
330 | |
---|
331 | ! for now, set excessoro=0.0 |
---|
332 | do jy=0,nyn(l)-1 |
---|
333 | do ix=0,nxn(l)-1 |
---|
334 | excessoron(ix,jy,l)=0.0 |
---|
335 | end do |
---|
336 | end do |
---|
337 | do jy=1,nyn(l)-2 |
---|
338 | do ix=1,nxn(l)-2 |
---|
339 | m=oron(ix,jy,l)+oron(ix-1,jy,l)+oron(ix+1,jy,l)+ & |
---|
340 | oron(ix,jy-1,l)+oron(ix,jy+1,l) |
---|
341 | m=m/5. |
---|
342 | excessoron(ix,jy,l)=(oron(ix,jy,l)-m)**2.+(oron(ix-1,jy,l)-m)**2. & |
---|
343 | +(oron(ix+1,jy,l)-m)**2.+(oron(ix,jy-1,l)-m)**2.+(oron(ix,jy+1,l)-m)**2. |
---|
344 | excessoron(ix,jy,l)=(excessoron(ix,jy,l)/5.)**0.5 |
---|
345 | end do |
---|
346 | end do |
---|
347 | |
---|
348 | |
---|
349 | ! |
---|
350 | ! identify the parent grid (which is probably "l-1", so this code |
---|
351 | ! may be more complicated that necessary) |
---|
352 | ! set xmet0n, ymet0n, which are the x,y coords of lower-left corner |
---|
353 | ! of a nested grid (in meters on the mother grid) |
---|
354 | ! |
---|
355 | ! note on dumc: |
---|
356 | ! the lower-left corner of the nested cell (0,0) coincides with the |
---|
357 | ! lower-left corner of the parent cell (i_parent_start,j_parent_start) |
---|
358 | ! this being the case, the center of nested cell (0,0) is shifted |
---|
359 | ! by (-dumc*parent_gridsize) relative to the center of the parent cell |
---|
360 | ! (for m_parent_grid_ratio = 2, 3, 4, 5; dumc = 1/4, 1/3, 3/8, 2/5) |
---|
361 | ! |
---|
362 | l_parent_nest_id(l) = -1 |
---|
363 | if (m_parent_grid_id(l) .gt. 0) then |
---|
364 | do lp = 0, l-1 |
---|
365 | if ( (l_parent_nest_id(l) .eq. -1) .and. & |
---|
366 | (m_parent_grid_id(l) .eq. m_grid_id(lp)) ) then |
---|
367 | |
---|
368 | l_parent_nest_id(l) = lp |
---|
369 | m = m_parent_grid_ratio(l) |
---|
370 | dumc = real(m-1)/real(m*2) |
---|
371 | if (lp .eq. 0) then |
---|
372 | xmet0n(l) = xmet0 + dx*(i_parent_start(l)-dumc) |
---|
373 | ymet0n(l) = ymet0 + dy*(j_parent_start(l)-dumc) |
---|
374 | else |
---|
375 | xmet0n(l) = xmet0n(lp) + dxn(lp)*(i_parent_start(l)-dumc) |
---|
376 | ymet0n(l) = ymet0n(lp) + dyn(lp)*(j_parent_start(l)-dumc) |
---|
377 | end if |
---|
378 | end if |
---|
379 | end do |
---|
380 | end if |
---|
381 | if (idiagbb .gt. 0) write(*,'(/a,3i8)') & |
---|
382 | 'l, m_grid_id(l), m_parent_grid_id(l)', & |
---|
383 | l, m_grid_id(l), m_parent_grid_id(l) |
---|
384 | |
---|
385 | if (l_parent_nest_id(l) .eq. -1) then |
---|
386 | write(*,'(/a,i3/)') & |
---|
387 | 'gridcheck_nests fatal error -- ' // & |
---|
388 | 'parent grid not found for l =', l |
---|
389 | stop |
---|
390 | end if |
---|
391 | |
---|
392 | ! |
---|
393 | ! diagnostics for testing the nesting calculations |
---|
394 | ! (set idiagbb=0 to turn it off) |
---|
395 | ! |
---|
396 | lp = l_parent_nest_id(l) |
---|
397 | if (idiagbb .gt. 0) then |
---|
398 | write(*,'(a,2i8)') 'l_parent, m_grid_id(l_parent) ', & |
---|
399 | lp, m_grid_id(lp) |
---|
400 | write(*,'(a,2i8)') 'm_parent_grid_ratio(l) ', & |
---|
401 | m_parent_grid_ratio(l) |
---|
402 | write(*,'(a,i8,f11.2)') & |
---|
403 | 'i_parent_start(l), xi_... ', & |
---|
404 | i_parent_start(l), i_parent_start(l)-dumc |
---|
405 | write(*,'(a,i8,f11.2)') & |
---|
406 | 'j_parent_start(l), yj_... ', & |
---|
407 | j_parent_start(l), j_parent_start(l)-dumc |
---|
408 | end if |
---|
409 | !23456789012345678901234567890123456789012345678901234567890123456789012 |
---|
410 | |
---|
411 | if (idiagbb .gt. 0) then |
---|
412 | write(*,*) |
---|
413 | do jy = j_parent_start(l)-1, j_parent_start(l) |
---|
414 | do ix = i_parent_start(l)-1, i_parent_start(l) |
---|
415 | if (lp .eq. 0) then |
---|
416 | write(*,'(a,2i7,2f11.4)') 'parent i,j,lon,lat', & |
---|
417 | ix, jy, xlon2d(ix,jy), ylat2d(ix,jy) |
---|
418 | else |
---|
419 | write(*,'(a,2i7,2f11.4)') 'parent i,j,lon,lat', & |
---|
420 | ix, jy, xlon2dn(ix,jy,lp), ylat2dn(ix,jy,lp) |
---|
421 | end if |
---|
422 | end do |
---|
423 | end do |
---|
424 | |
---|
425 | dumc = real( (m_parent_grid_ratio(l)-1) )/ & |
---|
426 | real( (m_parent_grid_ratio(l)*2) ) |
---|
427 | dumx = i_parent_start(l) - dumc |
---|
428 | dumy = j_parent_start(l) - dumc |
---|
429 | call xyindex_to_ll_wrf( lp, dumx, dumy, xaux1, yaux1 ) |
---|
430 | write(*,'(a,2f7.2,2f11.4)') 'par. xi,yj,lon,lat', & |
---|
431 | dumx, dumy, xaux1, yaux1 |
---|
432 | |
---|
433 | write(*,*) |
---|
434 | write(*,'(a,2i7,2f11.4)') 'nest i,j,lon,lat', & |
---|
435 | 0, 0, xlon2dn(0,0,l), ylat2dn(0,0,l) |
---|
436 | write(*,*) |
---|
437 | |
---|
438 | dumx = (xmet0n(l) - xmet0)/dx |
---|
439 | dumy = (ymet0n(l) - ymet0)/dy |
---|
440 | call xyindex_to_ll_wrf( 0, dumx, dumy, xaux1, yaux1 ) |
---|
441 | write(*,'(a,2f7.2,2f11.4)') 'mot. xi,yj,lon,lat', & |
---|
442 | dumx, dumy, xaux1, yaux1 |
---|
443 | |
---|
444 | iduma = max( 0, ifix(dumx) ) |
---|
445 | idumb = max( 0, ifix(dumy) ) |
---|
446 | do jy = idumb, idumb+1 |
---|
447 | do ix = iduma, iduma+1 |
---|
448 | write(*,'(a,2i7,2f11.4)') 'mother i,j,lon,lat', & |
---|
449 | ix, jy, xlon2d(ix,jy), ylat2d(ix,jy) |
---|
450 | end do |
---|
451 | end do |
---|
452 | end if |
---|
453 | |
---|
454 | |
---|
455 | |
---|
456 | ! Output of grid info |
---|
457 | !******************** |
---|
458 | |
---|
459 | write(*,'(/a,i2)') & |
---|
460 | 'gridcheck_nests -- nested domain #: ',l |
---|
461 | write(*,'(a,f12.1,a1,f12.1,a,f10.1)') & |
---|
462 | ' X coordinate range: ', & |
---|
463 | xmet0n(l),' to ',xmet0n(l)+(nxn(l)-1)*dxn(l), & |
---|
464 | ' Grid distance: ',dxn(l) |
---|
465 | write(*,'(a,f12.1,a1,f12.1,a,f10.1)') & |
---|
466 | ' Y coordinate range: ', & |
---|
467 | ymet0n(l),' to ',ymet0n(l)+(nyn(l)-1)*dyn(l), & |
---|
468 | ' Grid distance: ',dyn(l) |
---|
469 | write(*,*) |
---|
470 | |
---|
471 | |
---|
472 | ! Determine, how much the resolutions in the nests are enhanced as |
---|
473 | ! compared to the mother grid |
---|
474 | !***************************************************************** |
---|
475 | |
---|
476 | xresoln(l)=dx/dxn(l) |
---|
477 | yresoln(l)=dy/dyn(l) |
---|
478 | |
---|
479 | |
---|
480 | ! Determine the mother grid coordinates of the corner points of the |
---|
481 | ! nested grids |
---|
482 | ! Convert first to geographical coordinates, then to grid coordinates |
---|
483 | !******************************************************************** |
---|
484 | |
---|
485 | xaux1=xmet0n(l) |
---|
486 | xaux2=xmet0n(l)+real(nxn(l)-1)*dxn(l) |
---|
487 | yaux1=ymet0n(l) |
---|
488 | yaux2=ymet0n(l)+real(nyn(l)-1)*dyn(l) |
---|
489 | |
---|
490 | xln(l)=(xaux1-xmet0)/dx |
---|
491 | xrn(l)=(xaux2-xmet0)/dx |
---|
492 | yln(l)=(yaux1-ymet0)/dy |
---|
493 | yrn(l)=(yaux2-ymet0)/dy |
---|
494 | |
---|
495 | |
---|
496 | if ((xln(l).lt.0.).or.(yln(l).lt.0.).or. & |
---|
497 | (xrn(l).gt.real(nxmin1)).or.(yrn(l).gt.real(nymin1))) then |
---|
498 | write(*,*) 'gridcheck_nests error' |
---|
499 | write(*,*) 'Nested domain does not fit into mother domain' |
---|
500 | write(*,*) 'Execution is terminated.' |
---|
501 | stop |
---|
502 | endif |
---|
503 | |
---|
504 | |
---|
505 | ! check that the map projection routines are working |
---|
506 | call test_xyindex_to_ll_wrf( l ) |
---|
507 | |
---|
508 | |
---|
509 | ! Check, whether the heights of the model levels of the nested |
---|
510 | ! wind fields are consistent with those of the mother domain. |
---|
511 | ! If not, terminate model run. |
---|
512 | !************************************************************* |
---|
513 | |
---|
514 | ! first read eta_w_wrf_nest, eta_u_wrf_nest, and p_top_wrf_nest |
---|
515 | ! from the netcdf wrfout file |
---|
516 | itime = 1 |
---|
517 | |
---|
518 | varname = 'ZNW' |
---|
519 | do i = 1, ndims_max |
---|
520 | lendim_exp(i) = 0 |
---|
521 | lendim_max(i) = 1 |
---|
522 | end do |
---|
523 | lendim_exp(1) = nwz |
---|
524 | lendim_max(1) = nwzmax |
---|
525 | ndims_exp = 2 |
---|
526 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
527 | varname, eta_w_wrf_nest, & |
---|
528 | itime, & |
---|
529 | ndims, ndims_exp, ndims_max, & |
---|
530 | lendim, lendim_exp, lendim_max ) |
---|
531 | if (ierr .ne. 0) then |
---|
532 | write(*,*) |
---|
533 | write(*,*) '*** checkgrid -- error doing ncread of ZNW' |
---|
534 | stop |
---|
535 | end if |
---|
536 | |
---|
537 | varname = 'ZNU' |
---|
538 | do i = 1, ndims_max |
---|
539 | lendim_exp(i) = 0 |
---|
540 | lendim_max(i) = 1 |
---|
541 | end do |
---|
542 | lendim_exp(1) = nwz-1 |
---|
543 | lendim_max(1) = nwzmax |
---|
544 | ndims_exp = 2 |
---|
545 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
546 | varname, eta_u_wrf_nest, & |
---|
547 | itime, & |
---|
548 | ndims, ndims_exp, ndims_max, & |
---|
549 | lendim, lendim_exp, lendim_max ) |
---|
550 | if (ierr .ne. 0) then |
---|
551 | write(*,*) |
---|
552 | write(*,*) '*** checkgrid -- error doing ncread of ZNU' |
---|
553 | stop |
---|
554 | end if |
---|
555 | |
---|
556 | varname = 'P_TOP' |
---|
557 | do i = 1, ndims_max |
---|
558 | lendim_exp(i) = 0 |
---|
559 | lendim_max(i) = 1 |
---|
560 | end do |
---|
561 | lendim_exp(1) = 1 |
---|
562 | lendim_max(1) = 6 |
---|
563 | ndims_exp = 2 |
---|
564 | if (ext_scalar .lt. 0) ndims_exp=1 |
---|
565 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
566 | varname, p_top_wrf_nest, & |
---|
567 | itime, & |
---|
568 | ndims, ndims_exp, ndims_max, & |
---|
569 | lendim, lendim_exp, lendim_max ) |
---|
570 | if (ierr .ne. 0) then |
---|
571 | write(*,*) |
---|
572 | write(*,*) '*** checkgrid -- error doing ncread of P_TOP' |
---|
573 | stop |
---|
574 | end if |
---|
575 | |
---|
576 | |
---|
577 | do k = 1, nwz |
---|
578 | duma = 3.0e-7*max( abs(eta_w_wrf(k)), 1.0e-30 ) |
---|
579 | if (abs(eta_w_wrf(k)-eta_w_wrf_nest(k)) .gt. duma) then |
---|
580 | write(*,*) & |
---|
581 | 'FLEXPART gridcheck_nests error for nesting level',l |
---|
582 | write(*,*) & |
---|
583 | 'eta_w_wrf are not consistent with the mother domain' |
---|
584 | write(*,*) 'k, eta_w_wrf(k), eta_w_wrf_nest(k) =', & |
---|
585 | k, eta_w_wrf(k), eta_w_wrf_nest(k) |
---|
586 | stop |
---|
587 | endif |
---|
588 | end do |
---|
589 | |
---|
590 | do k = 1, nwz-1 |
---|
591 | duma = 3.0e-7*max( abs(eta_u_wrf(k)), 1.0e-30 ) |
---|
592 | if (abs(eta_u_wrf(k)-eta_u_wrf_nest(k)) .gt. duma) then |
---|
593 | write(*,*) & |
---|
594 | 'FLEXPART gridcheck_nests error for nesting level',l |
---|
595 | write(*,*) & |
---|
596 | 'eta_u_wrf are not consistent with the mother domain' |
---|
597 | write(*,*) 'k, eta_u_wrf(k), eta_u_wrf_nest(k) =', & |
---|
598 | k, eta_u_wrf(k), eta_u_wrf_nest(k) |
---|
599 | stop |
---|
600 | endif |
---|
601 | end do |
---|
602 | |
---|
603 | duma = 3.0e-7*max( abs(p_top_wrf), 1.0e-30 ) |
---|
604 | if (abs(p_top_wrf-p_top_wrf_nest) .gt. duma) then |
---|
605 | write(*,*) & |
---|
606 | 'FLEXPART gridcheck_nests error for nesting level',l |
---|
607 | write(*,*) & |
---|
608 | 'p_top_wrf are not consistent with the mother domain' |
---|
609 | write(*,*) 'p_top_wrf, p_top_wrf_nest', & |
---|
610 | p_top_wrf, p_top_wrf_nest |
---|
611 | stop |
---|
612 | endif |
---|
613 | |
---|
614 | |
---|
615 | ! done with nest l |
---|
616 | enddo |
---|
617 | |
---|
618 | return |
---|
619 | |
---|
620 | |
---|
621 | 999 write(*,*) |
---|
622 | write(*,*) ' ###########################################'// & |
---|
623 | '###### ' |
---|
624 | write(*,*) ' FLEXPART_WRF subroutine gridcheck_nests: ' // & |
---|
625 | 'nesting level ', l |
---|
626 | write(*,*) ' can not open input data file ' |
---|
627 | write(*,*) ' '//fnamenc |
---|
628 | write(*,*) ' or, an error occured in subr. read_ncwrfout_gridinfo' |
---|
629 | write(*,*) ' with ierr =', ierr |
---|
630 | write(*,*) ' ###########################################'// & |
---|
631 | '###### ' |
---|
632 | stop |
---|
633 | |
---|
634 | end subroutine gridcheck_nests |
---|
635 | |
---|
636 | |
---|