| 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 | |
|---|
