Ticket #14: gridcheck.f90

File gridcheck.f90, 22.6 KB (added by jbrioude, 11 years ago)

gridcheck version 3

Line 
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 = 100
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      lendim_exp(3) = 1
208      lendim_max(3) = 1
209      ndims_exp = 3
210      itime=1
211      print*, 'varname before',varname
212      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
213          varname, m_x(0,0,1), &
214          itime, &
215          ndims, ndims_exp, ndims_max, &
216          lendim, lendim_exp, lendim_max )
217      if (ierr .ne. 0) then
218      varname = 'MAPFAC_M'
219      lendim_exp(1) = nx
220      lendim_max(1) = nxmax
221      lendim_exp(2) = ny
222      lendim_max(2) = nymax
223
224      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
225          varname, m_x(0,0,1), &
226          itime, &
227          ndims, ndims_exp, ndims_max, &
228          lendim, lendim_exp, lendim_max )
229      endif
230      if (ierr .ne. 0) then
231          print*,'error doing MAP X'
232      varname = 'MAPFAC_UX'
233      lendim_exp(1) = nx+1
234      lendim_max(1) = nxmax
235      lendim_exp(2) = ny
236      lendim_max(2) = nymax
237      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
238          varname, m_u(0,0,1), &
239          itime, &
240          ndims, ndims_exp, ndims_max, &
241          lendim, lendim_exp, lendim_max )
242      do j = 0, nymin1
243      do i = 0, nxmin1
244      m_x(i,j,1)=(m_u(i,j,1)+m_u(i+1,j,1))*0.5
245      enddo
246      enddo
247      if (ierr .ne. 0) then
248          print*,'error doing MAP U'
249          print*,'NO MAP FACTOR IS GOING TO BE USED.'
250          print*,'LARGE UNCERTAINTIES TO BE EXPECTED'
251      do j = 0, nymin1
252      do i = 0, nxmin1
253      m_x(i,j,1)=1.
254      enddo
255      enddo
256      end if
257      end if
258!     do j = 0, nymin1
259!     do i = 0, nxmin1
260!     m_x(i,j,1)=1.
261!     enddo
262!     enddo
263
264      varname = 'MAPFAC_MY'
265      lendim_exp(1) = nx
266      lendim_max(1) = nxmax
267      lendim_exp(2) = ny
268      lendim_max(2) = nymax
269
270      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
271          varname, m_y(0,0,1), &
272          itime, &
273          ndims, ndims_exp, ndims_max, &
274          lendim, lendim_exp, lendim_max )
275      if (ierr .ne. 0) then
276      varname = 'MAPFAC_M'
277      lendim_exp(1) = nx
278      lendim_max(1) = nxmax
279      lendim_exp(2) = ny
280      lendim_max(2) = nymax
281
282      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
283          varname, m_y(0,0,1), &
284          itime, &
285          ndims, ndims_exp, ndims_max, &
286          lendim, lendim_exp, lendim_max )
287      endif
288      if (ierr .ne. 0) then
289          print*,'error doing MAP Y'
290      varname = 'MAPFAC_VY'
291      lendim_exp(1) = nx
292      lendim_max(1) = nxmax
293      lendim_exp(2) = ny+1
294      lendim_max(2) = nymax
295      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
296          varname, m_v(0,0,1), &
297          itime, &
298          ndims, ndims_exp, ndims_max, &
299          lendim, lendim_exp, lendim_max )
300      do j = 0, nymin1
301      do i = 0, nxmin1
302      m_y(i,j,1)=(m_v(i,j,1)+m_v(i,j+1,1))*0.5
303      enddo
304      enddo
305      if (ierr .ne. 0) then
306          print*,'ERROR doing MAP V'
307          print*,'NO MAP FACTOR IS GOING TO BE USED.'
308          print*,'LARGE UNCERTAINTIES TO BE EXPECTED'
309      do j = 0, nymin1
310      do i = 0, nxmin1
311      m_y(i,j,1)=1.
312      enddo
313      enddo
314      end if
315      end if
316      lendim_exp(1) = nx
317      lendim_max(1) = nxmax
318      lendim_exp(2) = ny
319      lendim_max(2) = nymax
320!     do j = 0, nymin1
321!     do i = 0, nxmin1
322!     m_y(i,j,1)=1.
323!     enddo
324!     enddo
325
326!
327      idiagaa = 100
328
329      varname = 'XLAT'
330      do i = 1, ndims_max
331          lendim_exp(i) = 0
332          lendim_max(i) = 1
333      end do
334      itime = 1
335      lendim_exp(1) = nx
336      lendim_max(1) = nxmax
337      lendim_exp(2) = ny
338      lendim_max(2) = nymax
339      ndims_exp = 3
340      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
341          varname, ylat2d, &
342          itime, &
343          ndims, ndims_exp, ndims_max, &
344          lendim, lendim_exp, lendim_max )
345      if (ierr .ne. 0) then
346          write(*,*)
347          write(*,*) '*** checkgrid -- error doing ncread of XLAT'
348          stop
349      end if
350!       print*,'values',ylat2d(309,101),ylat2d(300,101),ylat2d(250,101)
351!       print*,'values',ylat2d(309,101),ylat2d(300,101),ylat2d(250,101)
352      varname = 'XLONG'
353      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
354          varname, xlon2d, &
355          itime, &
356          ndims, ndims_exp, ndims_max, &
357          lendim, lendim_exp, lendim_max )
358      if (ierr .ne. 0) then
359          write(*,*)
360          write(*,*) '*** checkgrid -- error doing ncread of XLONG'
361          stop
362      end if
363
364      varname = 'HGT'
365      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
366          varname, oro, &
367          itime, &
368          ndims, ndims_exp, ndims_max, &
369          lendim, lendim_exp, lendim_max )
370      if (ierr .ne. 0) then
371          write(*,*)
372          write(*,*) '*** checkgrid -- error doing ncread of HGT'
373          stop
374      end if
375
376! lsm = landsea mask == land fraction (or non-ocean fraction)
377! for now, set lsm=1.0 which means land
378      do jy=0,ny-1
379      do ix=0,nxfield-1
380!         lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
381          lsm(ix,jy)=1.0
382      end do
383      end do
384
385! for now, set excessoro=0.0
386      do jy=0,ny-1
387      do ix=0,nxfield-1
388          excessoro(ix,jy)=0.0
389      end do
390      end do
391      do jy=1,ny-2
392      do ix=1,nxfield-2
393      m=4.*oro(ix,jy)+oro(ix-1,jy)+oro(ix+1,jy)+oro(ix,jy-1)+oro(ix,jy+1)
394      m=m/8.
395      excessoro(ix,jy)=4.*(oro(ix,jy)-m)**2.+(oro(ix-1,jy)-m)**2. &
396      +(oro(ix+1,jy)-m)**2.+(oro(ix,jy-1)-m)**2.+(oro(ix,jy+1)-m)**2.
397      excessoro(ix,jy)=(excessoro(ix,jy)/8.)**0.5
398      end do
399      end do
400
401
402! check that the map projection routines are working
403      call test_xyindex_to_ll_wrf( 0 )
404
405
406! If desired, shift all grids by nxshift grid cells
407!**************************************************
408
409!     if (xglobal) then
410!       call shift_field_0(oro,nxfield,ny)
411!       call shift_field_0(lsm,nxfield,ny)
412!       call shift_field_0(excessoro,nxfield,ny)
413!     endif
414
415
416! CALCULATE VERTICAL DISCRETIZATION OF WRF MODEL
417! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
418
419! read eta_w_wrf, eta_u_wrf, and p_top_wrf from the netcdf wrfout file
420      itime = 1
421
422      varname = 'ZNW'
423      do i = 1, ndims_max
424          lendim_exp(i) = 0
425          lendim_max(i) = 1
426      end do
427      lendim_exp(1) = nwz
428      lendim_max(1) = nwzmax
429      ndims_exp = 2
430      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
431          varname, eta_w_wrf, &
432          itime, &
433          ndims, ndims_exp, ndims_max, &
434          lendim, lendim_exp, lendim_max )
435      if (ierr .ne. 0) then
436        fnamenc2='wrfout_d03_zn.nc'
437      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, &
438          varname, eta_w_wrf, &
439          itime, &
440          ndims, ndims_exp, ndims_max, &
441          lendim, lendim_exp, lendim_max )
442      if (ierr .ne. 0) then
443
444          write(*,*)
445          write(*,*) '*** checkgrid -- error doing ncread of ZNW'
446!         stop
447      end if
448      end if
449
450      varname = 'ZNU'
451      do i = 1, ndims_max
452          lendim_exp(i) = 0
453          lendim_max(i) = 1
454      end do
455      lendim_exp(1) = nwz-1
456      lendim_max(1) = nwzmax
457      ndims_exp = 2
458      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
459          varname, eta_u_wrf, &
460          itime, &
461          ndims, ndims_exp, ndims_max, &
462          lendim, lendim_exp, lendim_max )
463      if (ierr .ne. 0) then
464        fnamenc2='wrfout_d03_zn.nc'
465      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc2, &
466          varname, eta_u_wrf, &
467          itime, &
468          ndims, ndims_exp, ndims_max, &
469          lendim, lendim_exp, lendim_max )
470       
471      if (ierr .ne. 0) then
472          write(*,*)
473          write(*,*) '*** checkgrid -- error doing ncread of ZNU'
474          stop
475      end if
476      end if
477
478      varname = 'P_TOP'
479      do i = 1, ndims_max
480          lendim_exp(i) = 0
481          lendim_max(i) = 1
482      end do
483      lendim_exp(1) = 1
484      lendim_max(1) = 6
485      ndims_exp = 2
486      if (ext_scalar .lt. 0) ndims_exp = 1
487      call read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
488          varname, p_top_wrf, &
489          itime, &
490          ndims, ndims_exp, ndims_max, &
491          lendim, lendim_exp, lendim_max )
492      if (ierr .ne. 0) then
493          write(*,*)
494          write(*,*) '*** checkgrid -- error doing ncread of P_TOP'
495          stop
496      end if
497
498! diagnostics for testing
499      if (idiagaa .gt. 0) then
500          write(*,*)
501          write(*,*) 'k, eta_w_wrf, eta_u_wrf ='
502          write(*,'(i3,2f11.6)')  &
503              (kz, eta_w_wrf(kz), eta_u_wrf(kz), kz=1,nwz-1)
504          kz = nwz
505          write(*,'(i3,2f11.6)') kz, eta_w_wrf(kz)
506          write(*,*)
507          write(*,*) 'p_top_wrf =', p_top_wrf
508          write(*,*)
509
510          duma = 0.0
511          dumb = 1.0e30
512          dumc = -1.0e30
513          do jy = 0, ny-1
514          do ix = 0, nx-1
515              duma = duma + oro(ix,jy)
516              dumb = min( dumb, oro(ix,jy) )
517              dumc = max( dumc, oro(ix,jy) )
518          end do
519          end do
520          duma = duma/(nx*ny)
521          write(*,*) 'oro avg, min, max =', duma, dumb, dumc
522          write(*,*)
523      end if
524
525
526!
527! the wrf eta vertical grid at layer boundaries (w grid) and
528! layer centers (u grid) is defined by
529!       eta_w_wrf(kz) = (pdh_w(kz) - p_top_wrf)/(pdh_surface - p_top_wrf)
530!       eta_u_wrf(kz) = (pdh_u(kz) - p_top_wrf)/(pdh_surface - p_top_wrf)
531! where "pdh_" refers to the dry hydrostatic component of the pressure
532!
533! so
534!       pdh_w(kz) = ((1.0 - eta_w_wrf(kz))*p_top_wrf) + eta_w_wrf(kz)*pdh_surface
535!       pdh_u(kz) = ((1.0 - eta_u_wrf(kz))*p_top_wrf) + eta_u_wrf(kz)*pdh_surface
536!
537! the ecmwf eta vertical grid is defined by
538!       p_w(kz) = akm(kz) + bkm(kz)*p_surface
539!       p_u(kz) = akz(kz) + bkz(kz)*p_surface
540!
541! the following definitions of akm, bkm, akz, bkz for wrf would be roughly
542! consistent those for ecmwf EXCEPT that for wrf, they involve the
543! dry hydrostatic component of the pressure
544!     do kz = 1, nwz
545!         akm(kz) = (1.0 - eta_w_wrf(kz))*p_top_wrf
546!         bkm(kz) = eta_w_wrf(kz)
547!     end do
548!     do kz = 1, nuvz
549!         akz(kz) = (1.0 - eta_u_wrf(kz))*p_top_wrf
550!         bkz(kz) = eta_u_wrf(kz)
551!     end do
552!
553! *** in FLEXPART_WRF we decided to used pressure from the met. files
554!     and drop the akz/bkz/akm/bkm entirely ***
555!
556
557
558! in FLEXPART_ECMWF, the U, V, & T-grid levels are always shifted up by 1,
559!    and an extra near-surface level is defined at kz=1
560!    which is loaded with the 10 m winds and 2 m temperature & humidity
561! for FLEXPART_WRF, this is optional, and is done when add_sfc_level=1
562! (Note -- it will take a lot of effort to get rid of this augmented
563!    level because many of the surface & boundary layer routines
564!    are expecting it.  so for now, always augment.)
565
566      dump1 = (101325.0-p_top_wrf)*eta_w_wrf(1) + p_top_wrf
567      dump2 = (101325.0-p_top_wrf)*eta_w_wrf(2) + p_top_wrf
568      dumdz = log(dump1/dump2)*8.4e3
569      write(*,*)
570!     write(*,*) 'add_sfc_level =', add_sfc_level
571!     write(*,*) 'WRF layer 1 approx. thickness =', dumdz
572
573      if (add_sfc_level .eq. 1) then
574          nuvz = nuvz + 1
575      else
576          write(*,'(/a/a/)') '*** gridcheck fatal error ***', &
577              '    add_sfc_level=0 is not yet implemented'
578!         stop
579      end if
580
581
582!*******************************************************************************
583! following comments are from FLEXPART_ECMWF.  This options for doubled vertical
584! resolution has not been tried in FLEXPART_WRF, but it probably could be done
585! with little effort.
586! ------------------------------------------------------------------------------
587! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
588! upon the transformation to z levels. In order to save computer memory, this is
589! not done anymore in the standard version. However, this option can still be
590! switched on by replacing the following lines with those below, that are
591! currently commented out. For this, similar changes are necessary in
592! verttransform.f and verttranform_nests.f
593!*******************************************************************************
594
595      nz=nuvz
596      if (nz.gt.nzmax) stop 'nzmax too small'
597
598!     do 100 i=1,nuvz
599!       aknew(i)=akz(i)
600!100    bknew(i)=bkz(i)
601
602! Switch on following lines to use doubled vertical resolution
603!*************************************************************
604!     nz=nuvz+nwz-1
605!     if (nz.gt.nzmax) stop 'nzmax too small'
606!     do 100 i=1,nwz
607!       aknew(2*(i-1)+1)=akm(i)
608!00     bknew(2*(i-1)+1)=bkm(i)
609!     do 110 i=2,nuvz
610!       aknew(2*(i-1))=akz(i)
611!10     bknew(2*(i-1))=bkz(i)
612! End doubled vertical resolution
613
614
615
616! Determine the uppermost level for which the convection scheme shall be applied
617! by assuming that there is no convection above 50 hPa (for standard SLP)
618!
619! FLEXPART_WRF - use approx. pressures to set nconvlev, and limit it to nuvz-2
620!*******************************************************************************
621
622      do i=1,nuvz-2
623!     do i=1,nuvz-1
624!       pint=akz(i)+bkz(i)*101325.
625        pint = (101325.0-p_top_wrf)*eta_u_wrf(i) + p_top_wrf
626        if (pint.lt.5000.0) goto 96
627      enddo
62896    nconvlev=i
629      nconvlev = min( nconvlev, nuvz-2 )
630      if (nconvlev.gt.nconvlevmax-1) then
631        nconvlev=nconvlevmax-1
632        pint = (101325.0-p_top_wrf)*eta_u_wrf(nconvlev) + p_top_wrf
633        write(*,*) 'Attention, convection only calculated up to ', &
634            pint*0.01, ' hPa'
635      endif
636
637
638! Output of grid info
639!********************
640
641      write(*,*)
642      write(*,*)
643      write(*,'(a/a,2i7/a,2i7//a,3i7/a,2i7/a,4i7)')  &
644        '# of vertical levels in WRF data',  &
645        '    n_bottom_top & "true" nuvz:', n_bottom_top,  &
646                                           (nuvz-add_sfc_level), &
647        '    nwz &     "augmented" nuvz:', nwz, nuvz, &
648        '    nwzmax, nuvzmax, nzmax    :', nwzmax, nuvzmax, nzmax, &
649        '    nconvlevmax, nconvlev     :', nconvlevmax, nconvlev, &
650        '    nx, ny, nxmax, nymax      :', nx, ny, nxmax, nymax 
651      write(*,*)
652      write(*,'(a)') 'Mother domain:'
653      write(*,'(a,f10.1,a1,f10.1,a,f10.1)') '  east-west   range: ', &
654        xmet0,' to ',xmet0+(nx-1)*dx,'   Grid distance: ',dx
655      write(*,'(a,f10.1,a1,f10.1,a,f10.1)') '  south-north range: ', &
656        ymet0,' to ',ymet0+(ny-1)*dy,'   Grid distance: ',dy
657      write(*,*)
658
659
660!
661! all done
662!
663      return
664
665
666! file open error
667999   write(*,*) 
668      write(*,*) ' ###########################################'// &
669                 '###### '
670      write(*,*) '       TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
671      write(*,*) ' CAN NOT OPEN INPUT DATA FILE = '
672      write(*,*) wfname(ifn)
673      write(*,*) ' ###########################################'// &
674                 '###### '
675      write(*,*)
676      stop
677
678end subroutine gridcheck
679
680
681
hosted by ZAMG