1 | subroutine gridcheck |
---|
2 | !********************************************************************** |
---|
3 | ! * |
---|
4 | ! TRAJECTORY MODEL SUBROUTINE GRIDCHECK * |
---|
5 | ! FLEXPART VERSION -> DO NOT USE IN FLEXTRA, PRAEPRO * |
---|
6 | ! * |
---|
7 | !********************************************************************** |
---|
8 | ! * |
---|
9 | ! AUTHOR: R. Easter & J. Fast, PNNL * |
---|
10 | ! DATE: 2005-autumn-?? * |
---|
11 | ! * |
---|
12 | ! Dec 2005, R. Easter - changed names of "*lon0*" & "*lat0*" variables* |
---|
13 | ! * |
---|
14 | !********************************************************************** |
---|
15 | ! * |
---|
16 | ! DESCRIPTION: * |
---|
17 | ! * |
---|
18 | ! Note: This is the FLEXPART_WRF version of subroutine gridcheck. * |
---|
19 | ! The computational grid is the WRF x-y grid rather than lat-lon. * |
---|
20 | ! There are many differences from the FLEXPART version. * |
---|
21 | ! * |
---|
22 | ! This subroutine determines the grid specifications * |
---|
23 | ! of the WRF model run from the first met. input file. * |
---|
24 | ! (longitudes & latitudes, number of grid points, grid distance, * |
---|
25 | ! vertical discretization) * |
---|
26 | ! The consistancy of met. input files is checked in the routine * |
---|
27 | ! "readwind" each call. (Grid info must not change.) * |
---|
28 | ! * |
---|
29 | ! * |
---|
30 | ! XMET0,YMET0 -- x,y coordinates (in m) of the lower left * |
---|
31 | ! "T-grid" point * |
---|
32 | ! NOTE: These replace XLON0,YLAT0 of FLEXPART (ECMWF), * |
---|
33 | ! which uses longitude,longitude (degrees) as coordinates. * |
---|
34 | ! * |
---|
35 | ! NX number of grid points x-direction * |
---|
36 | ! NY number of grid points y-direction * |
---|
37 | ! DX grid distance x-direction (in m) * |
---|
38 | ! DY grid distance y-direction (in m) * |
---|
39 | ! NUVZ number of grid points for horizontal wind * |
---|
40 | ! components in z direction * |
---|
41 | ! NWZ number of grid points for vertical wind * |
---|
42 | ! component in z direction * |
---|
43 | ! * |
---|
44 | ! For WRF, the "W" grid has NWZ == "bottom_top_stag" levels * |
---|
45 | ! For WRF, the "U", "V", and "T" grids have * |
---|
46 | ! NWZ-1 == "bottom_top" levels. * |
---|
47 | ! In the ecmwf FLEXPART, the "U", "V", and "T" grids are * |
---|
48 | ! "augmented" with an additional "surface" layer, * |
---|
49 | ! and thus have NUVZ==NWZ levels * |
---|
50 | ! Because of the high vertical resolution often used in WRF, * |
---|
51 | ! it may be desirable to eliminate this "surface layer". * |
---|
52 | ! * |
---|
53 | ! * |
---|
54 | ! UVHEIGHT(1)- heights of gridpoints where u and v are * |
---|
55 | ! UVHEIGHT(NUVZ) given * |
---|
56 | ! WHEIGHT(1)- heights of gridpoints where w is given * |
---|
57 | ! WHEIGHT(NWZ) * |
---|
58 | ! * |
---|
59 | !********************************************************************** |
---|
60 | ! |
---|
61 | ! include 'includepar' |
---|
62 | ! include 'includecom' |
---|
63 | ! include 'includeconv' |
---|
64 | ! use grib_api |
---|
65 | use par_mod |
---|
66 | use com_mod |
---|
67 | use conv_mod |
---|
68 | use cmapf_mod, only: stlmbr,stcm2p |
---|
69 | |
---|
70 | integer,parameter :: ndims_max=4 |
---|
71 | integer :: i, idiagaa, ierr, ifn, itime, ix |
---|
72 | integer :: jy |
---|
73 | integer :: kz |
---|
74 | integer :: lendim(ndims_max), lendim_exp(ndims_max), & |
---|
75 | lendim_max(ndims_max) |
---|
76 | integer :: ndims, ndims_exp, & |
---|
77 | ext_scalar,pbl_physics,mp_physics_dum |
---|
78 | integer :: n_west_east, n_south_north, n_bottom_top |
---|
79 | |
---|
80 | real :: dx_met, dy_met |
---|
81 | real :: duma, dumb, dumc |
---|
82 | real :: dump1, dump2, dumdz |
---|
83 | real :: pint,m |
---|
84 | |
---|
85 | character(len=160) :: fnamenc, varname,fnamenc2 |
---|
86 | |
---|
87 | real(kind=4) :: m_u(0:nxmax-1,0:nymax-1,1) |
---|
88 | real(kind=4) :: m_v(0:nxmax-1,0:nymax-1,1) |
---|
89 | |
---|
90 | |
---|
91 | |
---|
92 | ! |
---|
93 | ! get grid info from the wrf netcdf file |
---|
94 | ! |
---|
95 | ! write(*,'(//a)') 'gridcheck output' |
---|
96 | |
---|
97 | if(ideltas.gt.0) then |
---|
98 | ifn=1 |
---|
99 | else |
---|
100 | ifn=numbwf |
---|
101 | endif |
---|
102 | fnamenc = path(2)(1:length(2))//wfname(ifn) |
---|
103 | idiagaa = 0 |
---|
104 | |
---|
105 | call read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, & |
---|
106 | n_west_east, n_south_north, n_bottom_top, & |
---|
107 | dx_met, dy_met, & |
---|
108 | m_grid_id(0), m_parent_grid_id(0), m_parent_grid_ratio(0), & |
---|
109 | i_parent_start(0), j_parent_start(0), & |
---|
110 | map_proj_id, map_stdlon, map_truelat1, map_truelat2, & |
---|
111 | ext_scalar,pbl_physics,mp_physics_dum ) |
---|
112 | |
---|
113 | mp_physics=mp_physics_dum |
---|
114 | |
---|
115 | if (ierr.ne.0) goto 999 |
---|
116 | |
---|
117 | ! |
---|
118 | ! set grid dimension and size variables |
---|
119 | ! |
---|
120 | nx = n_west_east |
---|
121 | nxfield = nx |
---|
122 | ny = n_south_north |
---|
123 | |
---|
124 | nxmin1=nx-1 |
---|
125 | nymin1=ny-1 |
---|
126 | ! print*,'nxo',nxmin1,nymin1 |
---|
127 | nuvz = n_bottom_top |
---|
128 | nwz = n_bottom_top + 1 |
---|
129 | nlev_ec = n_bottom_top |
---|
130 | |
---|
131 | ! for FLEXPART_WRF, x & y coords are in meters, and |
---|
132 | ! we define lower-left corner of outermost (mother) grid == (0.0,0.0) |
---|
133 | xmet0 = 0.0 |
---|
134 | ymet0 = 0.0 |
---|
135 | |
---|
136 | dx = dx_met |
---|
137 | dy = dy_met |
---|
138 | |
---|
139 | dxconst = 1.0/dx |
---|
140 | dyconst = 1.0/dy |
---|
141 | |
---|
142 | l_parent_nest_id(0) = -1 |
---|
143 | |
---|
144 | xglobal=.false. |
---|
145 | if (nxshift.ne.0) stop & |
---|
146 | 'nxshift (par_mod) must be zero for non-global domain' |
---|
147 | |
---|
148 | sglobal=.false. |
---|
149 | switchsouthg=999999. |
---|
150 | |
---|
151 | nglobal=.false. |
---|
152 | switchnorthg=999999. |
---|
153 | |
---|
154 | if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative' |
---|
155 | if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large' |
---|
156 | |
---|
157 | ! |
---|
158 | ! check that grid dimensions are not too big |
---|
159 | ! |
---|
160 | ! FLEXPART_WRF 07-nov-2005 - require (nx+1 .le. nxmax) and (ny+1 .le. nymax) |
---|
161 | ! because u,v in met. files are on staggered grid |
---|
162 | ! if (nx.gt.nxmax) then |
---|
163 | if (nx+1 .gt. nxmax) then |
---|
164 | write(*,*) 'FLEXPART error: Too many grid points in x direction.' |
---|
165 | write(*,*) 'Reduce resolution of wind fields.' |
---|
166 | write(*,*) 'Or change parameter settings in file par_mod.' |
---|
167 | write(*,*) 'nx+1, nxmax =', nx,nxmax |
---|
168 | stop |
---|
169 | endif |
---|
170 | |
---|
171 | ! if (ny.gt.nymax) then |
---|
172 | if (ny+1 .gt. nymax) then |
---|
173 | write(*,*) 'FLEXPART error: Too many grid points in y direction.' |
---|
174 | write(*,*) 'Reduce resolution of wind fields.' |
---|
175 | write(*,*) 'Or change parameter settings in file par_mod.' |
---|
176 | write(*,*) 'ny+1, nymax =', ny,nymax |
---|
177 | stop |
---|
178 | endif |
---|
179 | |
---|
180 | if (nuvz+1.gt.nuvzmax) then |
---|
181 | write(*,*) 'FLEXPART error: Too many u,v grid points in z '// & |
---|
182 | 'direction.' |
---|
183 | write(*,*) 'Reduce resolution of wind fields.' |
---|
184 | write(*,*) 'Or change parameter settings in file par_mod.' |
---|
185 | write(*,*) nuvz+1,nuvzmax |
---|
186 | stop |
---|
187 | endif |
---|
188 | |
---|
189 | if (nwz.gt.nwzmax) then |
---|
190 | write(*,*) 'FLEXPART error: Too many w grid points in z '// & |
---|
191 | 'direction.' |
---|
192 | write(*,*) 'Reduce resolution of wind fields.' |
---|
193 | write(*,*) 'Or change parameter settings in file par_mod.' |
---|
194 | write(*,*) nwz,nwzmax |
---|
195 | stop |
---|
196 | endif |
---|
197 | |
---|
198 | ! |
---|
199 | ! read latitude and longitude |
---|
200 | ! read oro, lsm, and excessoro |
---|
201 | |
---|
202 | varname = 'MAPFAC_MX' |
---|
203 | lendim_exp(1) = nx |
---|
204 | lendim_max(1) = nxmax |
---|
205 | lendim_exp(2) = ny |
---|
206 | lendim_max(2) = nymax |
---|
207 | ndims_exp = 3 |
---|
208 | itime=1 |
---|
209 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
210 | varname, m_x(0,0,1), & |
---|
211 | itime, & |
---|
212 | ndims, ndims_exp, ndims_max, & |
---|
213 | lendim, lendim_exp, lendim_max ) |
---|
214 | if (ierr .ne. 0) then |
---|
215 | varname = 'MAPFAC_M' |
---|
216 | lendim_exp(1) = nx |
---|
217 | lendim_max(1) = nxmax |
---|
218 | lendim_exp(2) = ny |
---|
219 | lendim_max(2) = nymax |
---|
220 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
221 | varname, m_x(0,0,1), & |
---|
222 | itime, & |
---|
223 | ndims, ndims_exp, ndims_max, & |
---|
224 | lendim, lendim_exp, lendim_max ) |
---|
225 | endif |
---|
226 | if (ierr .ne. 0) then |
---|
227 | print*,'error doing MAP X' |
---|
228 | varname = 'MAPFAC_UX' |
---|
229 | lendim_exp(1) = nx+1 |
---|
230 | lendim_max(1) = nxmax |
---|
231 | lendim_exp(2) = ny |
---|
232 | lendim_max(2) = nymax |
---|
233 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
234 | varname, m_u(0,0,1), & |
---|
235 | itime, & |
---|
236 | ndims, ndims_exp, ndims_max, & |
---|
237 | lendim, lendim_exp, lendim_max ) |
---|
238 | do j = 0, nymin1 |
---|
239 | do i = 0, nxmin1 |
---|
240 | m_x(i,j,1)=(m_u(i,j,1)+m_u(i+1,j,1))*0.5 |
---|
241 | enddo |
---|
242 | enddo |
---|
243 | if (ierr .ne. 0) then |
---|
244 | print*,'error doing MAP U' |
---|
245 | print*,'NO MAP FACTOR IS GOING TO BE USED.' |
---|
246 | print*,'LARGE UNCERTAINTIES TO BE EXPECTED' |
---|
247 | do j = 0, nymin1 |
---|
248 | do i = 0, nxmin1 |
---|
249 | m_x(i,j,1)=1. |
---|
250 | enddo |
---|
251 | enddo |
---|
252 | end if |
---|
253 | end if |
---|
254 | ! do j = 0, nymin1 |
---|
255 | ! do i = 0, nxmin1 |
---|
256 | ! m_x(i,j,1)=1. |
---|
257 | ! enddo |
---|
258 | ! enddo |
---|
259 | |
---|
260 | varname = 'MAPFAC_MY' |
---|
261 | lendim_exp(1) = nx |
---|
262 | lendim_max(1) = nxmax |
---|
263 | lendim_exp(2) = ny |
---|
264 | lendim_max(2) = nymax |
---|
265 | |
---|
266 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
267 | varname, m_y(0,0,1), & |
---|
268 | itime, & |
---|
269 | ndims, ndims_exp, ndims_max, & |
---|
270 | lendim, lendim_exp, lendim_max ) |
---|
271 | if (ierr .ne. 0) then |
---|
272 | varname = 'MAPFAC_M' |
---|
273 | lendim_exp(1) = nx |
---|
274 | lendim_max(1) = nxmax |
---|
275 | lendim_exp(2) = ny |
---|
276 | lendim_max(2) = nymax |
---|
277 | |
---|
278 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
279 | varname, m_y(0,0,1), & |
---|
280 | itime, & |
---|
281 | ndims, ndims_exp, ndims_max, & |
---|
282 | lendim, lendim_exp, lendim_max ) |
---|
283 | endif |
---|
284 | if (ierr .ne. 0) then |
---|
285 | print*,'error doing MAP Y' |
---|
286 | varname = 'MAPFAC_VY' |
---|
287 | lendim_exp(1) = nx |
---|
288 | lendim_max(1) = nxmax |
---|
289 | lendim_exp(2) = ny+1 |
---|
290 | lendim_max(2) = nymax |
---|
291 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
292 | varname, m_v(0,0,1), & |
---|
293 | itime, & |
---|
294 | ndims, ndims_exp, ndims_max, & |
---|
295 | lendim, lendim_exp, lendim_max ) |
---|
296 | do j = 0, nymin1 |
---|
297 | do i = 0, nxmin1 |
---|
298 | m_y(i,j,1)=(m_v(i,j,1)+m_v(i,j+1,1))*0.5 |
---|
299 | enddo |
---|
300 | enddo |
---|
301 | if (ierr .ne. 0) then |
---|
302 | print*,'ERROR doing MAP V' |
---|
303 | print*,'NO MAP FACTOR IS GOING TO BE USED.' |
---|
304 | print*,'LARGE UNCERTAINTIES TO BE EXPECTED' |
---|
305 | do j = 0, nymin1 |
---|
306 | do i = 0, nxmin1 |
---|
307 | m_y(i,j,1)=1. |
---|
308 | enddo |
---|
309 | enddo |
---|
310 | end if |
---|
311 | end if |
---|
312 | lendim_exp(1) = nx |
---|
313 | lendim_max(1) = nxmax |
---|
314 | lendim_exp(2) = ny |
---|
315 | lendim_max(2) = nymax |
---|
316 | ! do j = 0, nymin1 |
---|
317 | ! do i = 0, nxmin1 |
---|
318 | ! m_y(i,j,1)=1. |
---|
319 | ! enddo |
---|
320 | ! enddo |
---|
321 | |
---|
322 | ! |
---|
323 | idiagaa = 0 |
---|
324 | |
---|
325 | varname = 'XLAT' |
---|
326 | do i = 1, ndims_max |
---|
327 | lendim_exp(i) = 0 |
---|
328 | lendim_max(i) = 1 |
---|
329 | end do |
---|
330 | itime = 1 |
---|
331 | lendim_exp(1) = nx |
---|
332 | lendim_max(1) = nxmax |
---|
333 | lendim_exp(2) = ny |
---|
334 | lendim_max(2) = nymax |
---|
335 | ndims_exp = 3 |
---|
336 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
337 | varname, ylat2d, & |
---|
338 | itime, & |
---|
339 | ndims, ndims_exp, ndims_max, & |
---|
340 | lendim, lendim_exp, lendim_max ) |
---|
341 | if (ierr .ne. 0) then |
---|
342 | write(*,*) |
---|
343 | write(*,*) '*** checkgrid -- error doing ncread of XLAT' |
---|
344 | stop |
---|
345 | end if |
---|
346 | ! print*,'values',ylat2d(309,101),ylat2d(300,101),ylat2d(250,101) |
---|
347 | ! print*,'values',ylat2d(309,101),ylat2d(300,101),ylat2d(250,101) |
---|
348 | varname = 'XLONG' |
---|
349 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
350 | varname, xlon2d, & |
---|
351 | itime, & |
---|
352 | ndims, ndims_exp, ndims_max, & |
---|
353 | lendim, lendim_exp, lendim_max ) |
---|
354 | if (ierr .ne. 0) then |
---|
355 | write(*,*) |
---|
356 | write(*,*) '*** checkgrid -- error doing ncread of XLONG' |
---|
357 | stop |
---|
358 | end if |
---|
359 | |
---|
360 | varname = 'HGT' |
---|
361 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
362 | varname, oro, & |
---|
363 | itime, & |
---|
364 | ndims, ndims_exp, ndims_max, & |
---|
365 | lendim, lendim_exp, lendim_max ) |
---|
366 | if (ierr .ne. 0) then |
---|
367 | write(*,*) |
---|
368 | write(*,*) '*** checkgrid -- error doing ncread of HGT' |
---|
369 | stop |
---|
370 | end if |
---|
371 | |
---|
372 | ! lsm = landsea mask == land fraction (or non-ocean fraction) |
---|
373 | ! for now, set lsm=1.0 which means land |
---|
374 | do jy=0,ny-1 |
---|
375 | do ix=0,nxfield-1 |
---|
376 | ! lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1) |
---|
377 | lsm(ix,jy)=1.0 |
---|
378 | end do |
---|
379 | end do |
---|
380 | |
---|
381 | ! for now, set excessoro=0.0 |
---|
382 | do jy=0,ny-1 |
---|
383 | do ix=0,nxfield-1 |
---|
384 | excessoro(ix,jy)=0.0 |
---|
385 | end do |
---|
386 | end do |
---|
387 | do jy=1,ny-2 |
---|
388 | do ix=1,nxfield-2 |
---|
389 | m=4.*oro(ix,jy)+oro(ix-1,jy)+oro(ix+1,jy)+oro(ix,jy-1)+oro(ix,jy+1) |
---|
390 | m=m/8. |
---|
391 | excessoro(ix,jy)=4.*(oro(ix,jy)-m)**2.+(oro(ix-1,jy)-m)**2. & |
---|
392 | +(oro(ix+1,jy)-m)**2.+(oro(ix,jy-1)-m)**2.+(oro(ix,jy+1)-m)**2. |
---|
393 | excessoro(ix,jy)=(excessoro(ix,jy)/8.)**0.5 |
---|
394 | end do |
---|
395 | end do |
---|
396 | |
---|
397 | |
---|
398 | ! check that the map projection routines are working |
---|
399 | call test_xyindex_to_ll_wrf( 0 ) |
---|
400 | |
---|
401 | |
---|
402 | ! If desired, shift all grids by nxshift grid cells |
---|
403 | !************************************************** |
---|
404 | |
---|
405 | ! if (xglobal) then |
---|
406 | ! call shift_field_0(oro,nxfield,ny) |
---|
407 | ! call shift_field_0(lsm,nxfield,ny) |
---|
408 | ! call shift_field_0(excessoro,nxfield,ny) |
---|
409 | ! endif |
---|
410 | |
---|
411 | |
---|
412 | ! CALCULATE VERTICAL DISCRETIZATION OF WRF MODEL |
---|
413 | ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM |
---|
414 | |
---|
415 | ! read eta_w_wrf, eta_u_wrf, and p_top_wrf from the netcdf wrfout file |
---|
416 | itime = 1 |
---|
417 | |
---|
418 | varname = 'ZNW' |
---|
419 | do i = 1, ndims_max |
---|
420 | lendim_exp(i) = 0 |
---|
421 | lendim_max(i) = 1 |
---|
422 | end do |
---|
423 | lendim_exp(1) = nwz |
---|
424 | lendim_max(1) = nwzmax |
---|
425 | ndims_exp = 2 |
---|
426 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
427 | varname, eta_w_wrf, & |
---|
428 | itime, & |
---|
429 | ndims, ndims_exp, ndims_max, & |
---|
430 | lendim, lendim_exp, lendim_max ) |
---|
431 | if (ierr .ne. 0) then |
---|
432 | fnamenc2='wrfout_d03_zn.nc' |
---|
433 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, & |
---|
434 | varname, eta_w_wrf, & |
---|
435 | itime, & |
---|
436 | ndims, ndims_exp, ndims_max, & |
---|
437 | lendim, lendim_exp, lendim_max ) |
---|
438 | if (ierr .ne. 0) then |
---|
439 | |
---|
440 | write(*,*) |
---|
441 | write(*,*) '*** checkgrid -- error doing ncread of ZNW' |
---|
442 | ! stop |
---|
443 | end if |
---|
444 | end if |
---|
445 | |
---|
446 | varname = 'ZNU' |
---|
447 | do i = 1, ndims_max |
---|
448 | lendim_exp(i) = 0 |
---|
449 | lendim_max(i) = 1 |
---|
450 | end do |
---|
451 | lendim_exp(1) = nwz-1 |
---|
452 | lendim_max(1) = nwzmax |
---|
453 | ndims_exp = 2 |
---|
454 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
455 | varname, eta_u_wrf, & |
---|
456 | itime, & |
---|
457 | ndims, ndims_exp, ndims_max, & |
---|
458 | lendim, lendim_exp, lendim_max ) |
---|
459 | if (ierr .ne. 0) then |
---|
460 | fnamenc2='wrfout_d03_zn.nc' |
---|
461 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, & |
---|
462 | varname, eta_u_wrf, & |
---|
463 | itime, & |
---|
464 | ndims, ndims_exp, ndims_max, & |
---|
465 | lendim, lendim_exp, lendim_max ) |
---|
466 | |
---|
467 | if (ierr .ne. 0) then |
---|
468 | write(*,*) |
---|
469 | write(*,*) '*** checkgrid -- error doing ncread of ZNU' |
---|
470 | stop |
---|
471 | end if |
---|
472 | end if |
---|
473 | |
---|
474 | varname = 'P_TOP' |
---|
475 | do i = 1, ndims_max |
---|
476 | lendim_exp(i) = 0 |
---|
477 | lendim_max(i) = 1 |
---|
478 | end do |
---|
479 | lendim_exp(1) = 1 |
---|
480 | lendim_max(1) = 6 |
---|
481 | ndims_exp = 2 |
---|
482 | if (ext_scalar .lt. 0) ndims_exp = 1 |
---|
483 | call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, & |
---|
484 | varname, p_top_wrf, & |
---|
485 | itime, & |
---|
486 | ndims, ndims_exp, ndims_max, & |
---|
487 | lendim, lendim_exp, lendim_max ) |
---|
488 | if (ierr .ne. 0) then |
---|
489 | write(*,*) |
---|
490 | write(*,*) '*** checkgrid -- error doing ncread of P_TOP' |
---|
491 | stop |
---|
492 | end if |
---|
493 | |
---|
494 | ! diagnostics for testing |
---|
495 | if (idiagaa .gt. 0) then |
---|
496 | write(*,*) |
---|
497 | write(*,*) 'k, eta_w_wrf, eta_u_wrf =' |
---|
498 | write(*,'(i3,2f11.6)') & |
---|
499 | (kz, eta_w_wrf(kz), eta_u_wrf(kz), kz=1,nwz-1) |
---|
500 | kz = nwz |
---|
501 | write(*,'(i3,2f11.6)') kz, eta_w_wrf(kz) |
---|
502 | write(*,*) |
---|
503 | write(*,*) 'p_top_wrf =', p_top_wrf |
---|
504 | write(*,*) |
---|
505 | |
---|
506 | duma = 0.0 |
---|
507 | dumb = 1.0e30 |
---|
508 | dumc = -1.0e30 |
---|
509 | do jy = 0, ny-1 |
---|
510 | do ix = 0, nx-1 |
---|
511 | duma = duma + oro(ix,jy) |
---|
512 | dumb = min( dumb, oro(ix,jy) ) |
---|
513 | dumc = max( dumc, oro(ix,jy) ) |
---|
514 | end do |
---|
515 | end do |
---|
516 | duma = duma/(nx*ny) |
---|
517 | write(*,*) 'oro avg, min, max =', duma, dumb, dumc |
---|
518 | write(*,*) |
---|
519 | end if |
---|
520 | |
---|
521 | |
---|
522 | ! |
---|
523 | ! the wrf eta vertical grid at layer boundaries (w grid) and |
---|
524 | ! layer centers (u grid) is defined by |
---|
525 | ! eta_w_wrf(kz) = (pdh_w(kz) - p_top_wrf)/(pdh_surface - p_top_wrf) |
---|
526 | ! eta_u_wrf(kz) = (pdh_u(kz) - p_top_wrf)/(pdh_surface - p_top_wrf) |
---|
527 | ! where "pdh_" refers to the dry hydrostatic component of the pressure |
---|
528 | ! |
---|
529 | ! so |
---|
530 | ! pdh_w(kz) = ((1.0 - eta_w_wrf(kz))*p_top_wrf) + eta_w_wrf(kz)*pdh_surface |
---|
531 | ! pdh_u(kz) = ((1.0 - eta_u_wrf(kz))*p_top_wrf) + eta_u_wrf(kz)*pdh_surface |
---|
532 | ! |
---|
533 | ! the ecmwf eta vertical grid is defined by |
---|
534 | ! p_w(kz) = akm(kz) + bkm(kz)*p_surface |
---|
535 | ! p_u(kz) = akz(kz) + bkz(kz)*p_surface |
---|
536 | ! |
---|
537 | ! the following definitions of akm, bkm, akz, bkz for wrf would be roughly |
---|
538 | ! consistent those for ecmwf EXCEPT that for wrf, they involve the |
---|
539 | ! dry hydrostatic component of the pressure |
---|
540 | ! do kz = 1, nwz |
---|
541 | ! akm(kz) = (1.0 - eta_w_wrf(kz))*p_top_wrf |
---|
542 | ! bkm(kz) = eta_w_wrf(kz) |
---|
543 | ! end do |
---|
544 | ! do kz = 1, nuvz |
---|
545 | ! akz(kz) = (1.0 - eta_u_wrf(kz))*p_top_wrf |
---|
546 | ! bkz(kz) = eta_u_wrf(kz) |
---|
547 | ! end do |
---|
548 | ! |
---|
549 | ! *** in FLEXPART_WRF we decided to used pressure from the met. files |
---|
550 | ! and drop the akz/bkz/akm/bkm entirely *** |
---|
551 | ! |
---|
552 | |
---|
553 | |
---|
554 | ! in FLEXPART_ECMWF, the U, V, & T-grid levels are always shifted up by 1, |
---|
555 | ! and an extra near-surface level is defined at kz=1 |
---|
556 | ! which is loaded with the 10 m winds and 2 m temperature & humidity |
---|
557 | ! for FLEXPART_WRF, this is optional, and is done when add_sfc_level=1 |
---|
558 | ! (Note -- it will take a lot of effort to get rid of this augmented |
---|
559 | ! level because many of the surface & boundary layer routines |
---|
560 | ! are expecting it. so for now, always augment.) |
---|
561 | |
---|
562 | dump1 = (101325.0-p_top_wrf)*eta_w_wrf(1) + p_top_wrf |
---|
563 | dump2 = (101325.0-p_top_wrf)*eta_w_wrf(2) + p_top_wrf |
---|
564 | dumdz = log(dump1/dump2)*8.4e3 |
---|
565 | write(*,*) |
---|
566 | ! write(*,*) 'add_sfc_level =', add_sfc_level |
---|
567 | ! write(*,*) 'WRF layer 1 approx. thickness =', dumdz |
---|
568 | |
---|
569 | if (add_sfc_level .eq. 1) then |
---|
570 | nuvz = nuvz + 1 |
---|
571 | else |
---|
572 | write(*,'(/a/a/)') '*** gridcheck fatal error ***', & |
---|
573 | ' add_sfc_level=0 is not yet implemented' |
---|
574 | ! stop |
---|
575 | end if |
---|
576 | |
---|
577 | |
---|
578 | !******************************************************************************* |
---|
579 | ! following comments are from FLEXPART_ECMWF. This options for doubled vertical |
---|
580 | ! resolution has not been tried in FLEXPART_WRF, but it probably could be done |
---|
581 | ! with little effort. |
---|
582 | ! ------------------------------------------------------------------------------ |
---|
583 | ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled |
---|
584 | ! upon the transformation to z levels. In order to save computer memory, this is |
---|
585 | ! not done anymore in the standard version. However, this option can still be |
---|
586 | ! switched on by replacing the following lines with those below, that are |
---|
587 | ! currently commented out. For this, similar changes are necessary in |
---|
588 | ! verttransform.f and verttranform_nests.f |
---|
589 | !******************************************************************************* |
---|
590 | |
---|
591 | nz=nuvz |
---|
592 | if (nz.gt.nzmax) stop 'nzmax too small' |
---|
593 | |
---|
594 | ! do 100 i=1,nuvz |
---|
595 | ! aknew(i)=akz(i) |
---|
596 | !100 bknew(i)=bkz(i) |
---|
597 | |
---|
598 | ! Switch on following lines to use doubled vertical resolution |
---|
599 | !************************************************************* |
---|
600 | ! nz=nuvz+nwz-1 |
---|
601 | ! if (nz.gt.nzmax) stop 'nzmax too small' |
---|
602 | ! do 100 i=1,nwz |
---|
603 | ! aknew(2*(i-1)+1)=akm(i) |
---|
604 | !00 bknew(2*(i-1)+1)=bkm(i) |
---|
605 | ! do 110 i=2,nuvz |
---|
606 | ! aknew(2*(i-1))=akz(i) |
---|
607 | !10 bknew(2*(i-1))=bkz(i) |
---|
608 | ! End doubled vertical resolution |
---|
609 | |
---|
610 | |
---|
611 | |
---|
612 | ! Determine the uppermost level for which the convection scheme shall be applied |
---|
613 | ! by assuming that there is no convection above 50 hPa (for standard SLP) |
---|
614 | ! |
---|
615 | ! FLEXPART_WRF - use approx. pressures to set nconvlev, and limit it to nuvz-2 |
---|
616 | !******************************************************************************* |
---|
617 | |
---|
618 | do i=1,nuvz-2 |
---|
619 | ! do i=1,nuvz-1 |
---|
620 | ! pint=akz(i)+bkz(i)*101325. |
---|
621 | pint = (101325.0-p_top_wrf)*eta_u_wrf(i) + p_top_wrf |
---|
622 | if (pint.lt.5000.0) goto 96 |
---|
623 | enddo |
---|
624 | 96 nconvlev=i |
---|
625 | nconvlev = min( nconvlev, nuvz-2 ) |
---|
626 | if (nconvlev.gt.nconvlevmax-1) then |
---|
627 | nconvlev=nconvlevmax-1 |
---|
628 | pint = (101325.0-p_top_wrf)*eta_u_wrf(nconvlev) + p_top_wrf |
---|
629 | write(*,*) 'Attention, convection only calculated up to ', & |
---|
630 | pint*0.01, ' hPa' |
---|
631 | endif |
---|
632 | |
---|
633 | |
---|
634 | ! Output of grid info |
---|
635 | !******************** |
---|
636 | |
---|
637 | write(*,*) |
---|
638 | write(*,*) |
---|
639 | write(*,'(a/a,2i7/a,2i7//a,3i7/a,2i7/a,4i7)') & |
---|
640 | '# of vertical levels in WRF data', & |
---|
641 | ' n_bottom_top & "true" nuvz:', n_bottom_top, & |
---|
642 | (nuvz-add_sfc_level), & |
---|
643 | ' nwz & "augmented" nuvz:', nwz, nuvz, & |
---|
644 | ' nwzmax, nuvzmax, nzmax :', nwzmax, nuvzmax, nzmax, & |
---|
645 | ' nconvlevmax, nconvlev :', nconvlevmax, nconvlev, & |
---|
646 | ' nx, ny, nxmax, nymax :', nx, ny, nxmax, nymax |
---|
647 | write(*,*) |
---|
648 | write(*,'(a)') 'Mother domain:' |
---|
649 | write(*,'(a,f10.1,a1,f10.1,a,f10.1)') ' east-west range: ', & |
---|
650 | xmet0,' to ',xmet0+(nx-1)*dx,' Grid distance: ',dx |
---|
651 | write(*,'(a,f10.1,a1,f10.1,a,f10.1)') ' south-north range: ', & |
---|
652 | ymet0,' to ',ymet0+(ny-1)*dy,' Grid distance: ',dy |
---|
653 | write(*,*) |
---|
654 | |
---|
655 | |
---|
656 | ! |
---|
657 | ! all done |
---|
658 | ! |
---|
659 | return |
---|
660 | |
---|
661 | |
---|
662 | ! file open error |
---|
663 | 999 write(*,*) |
---|
664 | write(*,*) ' ###########################################'// & |
---|
665 | '###### ' |
---|
666 | write(*,*) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:' |
---|
667 | write(*,*) ' CAN NOT OPEN INPUT DATA FILE = ' |
---|
668 | write(*,*) wfname(ifn) |
---|
669 | write(*,*) ' ###########################################'// & |
---|
670 | '###### ' |
---|
671 | write(*,*) |
---|
672 | stop |
---|
673 | |
---|
674 | end subroutine gridcheck |
---|
675 | |
---|
676 | |
---|
677 | |
---|