| 1 | program plotgrids |
|---|
| 2 | |
|---|
| 3 | use map_utils |
|---|
| 4 | |
|---|
| 5 | implicit none |
|---|
| 6 | |
|---|
| 7 | ! Parameters |
|---|
| 8 | integer, parameter :: MAX_DOMAINS = 21 |
|---|
| 9 | |
|---|
| 10 | ! Variables |
|---|
| 11 | integer :: iproj_type, n_domains, io_form_output, dyn_opt |
|---|
| 12 | integer :: i, j, max_dom, funit, io_form_geogrid |
|---|
| 13 | integer :: interval_seconds |
|---|
| 14 | |
|---|
| 15 | integer, dimension(MAX_DOMAINS) :: parent_grid_ratio, parent_id, ixdim, jydim |
|---|
| 16 | integer, dimension(MAX_DOMAINS) :: i_parent_start, j_parent_start, & |
|---|
| 17 | s_we, e_we, s_sn, e_sn, & |
|---|
| 18 | start_year, start_month, start_day, start_hour, & |
|---|
| 19 | end_year, end_month, end_day, end_hour |
|---|
| 20 | |
|---|
| 21 | real :: known_lat, known_lon, stand_lon, truelat1, truelat2, known_x, known_y, & |
|---|
| 22 | dxkm, dykm, phi, lambda, ref_lat, ref_lon, ref_x, ref_y |
|---|
| 23 | real :: dx, dy |
|---|
| 24 | real :: ri, rj, rlats, rlons, rlate, rlone |
|---|
| 25 | real :: polat , rot |
|---|
| 26 | real :: rparent_gridpts |
|---|
| 27 | real :: xa,xb,ya,yb,xxa,xxy,yya,yyb |
|---|
| 28 | real :: xs, xe, ys, ye |
|---|
| 29 | integer :: jproj, jgrid, jlts, iusout, idot, ier |
|---|
| 30 | integer :: ltype , idom |
|---|
| 31 | |
|---|
| 32 | real, dimension(MAX_DOMAINS) :: parent_ll_x, parent_ll_y, parent_ur_x, parent_ur_y |
|---|
| 33 | |
|---|
| 34 | character (len=128) :: geog_data_path, opt_output_from_geogrid_path, opt_geogrid_tbl_path |
|---|
| 35 | character (len=128), dimension(MAX_DOMAINS) :: geog_data_res |
|---|
| 36 | character (len=128) :: map_proj |
|---|
| 37 | character (len=128), dimension(MAX_DOMAINS) :: start_date, end_date |
|---|
| 38 | character (len=3) :: wrf_core |
|---|
| 39 | character (len=1) :: gridtype |
|---|
| 40 | |
|---|
| 41 | logical :: do_tiled_output |
|---|
| 42 | integer :: debug_level |
|---|
| 43 | logical :: is_used |
|---|
| 44 | |
|---|
| 45 | type (proj_info) :: map_projection |
|---|
| 46 | |
|---|
| 47 | namelist /share/ wrf_core, max_dom, start_date, end_date, & |
|---|
| 48 | start_year, end_year, start_month, end_month, & |
|---|
| 49 | start_day, end_day, start_hour, end_hour, & |
|---|
| 50 | interval_seconds, & |
|---|
| 51 | io_form_geogrid, opt_output_from_geogrid_path, debug_level |
|---|
| 52 | namelist /geogrid/ parent_id, parent_grid_ratio, & |
|---|
| 53 | i_parent_start, j_parent_start, s_we, e_we, s_sn, e_sn, & |
|---|
| 54 | map_proj, ref_x, ref_y, ref_lat, ref_lon, & |
|---|
| 55 | truelat1, truelat2, stand_lon, dx, dy, & |
|---|
| 56 | geog_data_res, geog_data_path, opt_geogrid_tbl_path |
|---|
| 57 | |
|---|
| 58 | ! Set defaults for namelist variables |
|---|
| 59 | debug_level = 0 |
|---|
| 60 | io_form_geogrid = 2 |
|---|
| 61 | wrf_core = 'ARW' |
|---|
| 62 | max_dom = 1 |
|---|
| 63 | geog_data_path = 'NOT_SPECIFIED' |
|---|
| 64 | ref_x = NAN |
|---|
| 65 | ref_y = NAN |
|---|
| 66 | ref_lat = NAN |
|---|
| 67 | ref_lon = NAN |
|---|
| 68 | dx = 10000. |
|---|
| 69 | dy = 10000. |
|---|
| 70 | map_proj = 'Lambert' |
|---|
| 71 | truelat1 = NAN |
|---|
| 72 | truelat2 = NAN |
|---|
| 73 | stand_lon = NAN |
|---|
| 74 | do i=1,MAX_DOMAINS |
|---|
| 75 | geog_data_res(i) = 'default' |
|---|
| 76 | parent_id(i) = 1 |
|---|
| 77 | parent_grid_ratio(i) = INVALID |
|---|
| 78 | s_we(i) = 1 |
|---|
| 79 | e_we(i) = INVALID |
|---|
| 80 | s_sn(i) = 1 |
|---|
| 81 | e_sn(i) = INVALID |
|---|
| 82 | start_year(i) = 0 |
|---|
| 83 | start_month(i) = 0 |
|---|
| 84 | start_day(i) = 0 |
|---|
| 85 | start_hour(i) = 0 |
|---|
| 86 | end_year(i) = 0 |
|---|
| 87 | end_month(i) = 0 |
|---|
| 88 | end_day(i) = 0 |
|---|
| 89 | end_hour(i) = 0 |
|---|
| 90 | start_date(i) = '0000-00-00_00:00:00' |
|---|
| 91 | end_date(i) = '0000-00-00_00:00:00' |
|---|
| 92 | end do |
|---|
| 93 | opt_output_from_geogrid_path = './' |
|---|
| 94 | opt_geogrid_tbl_path = 'geogrid/' |
|---|
| 95 | interval_seconds = INVALID |
|---|
| 96 | |
|---|
| 97 | ! Read parameters from Fortran namelist |
|---|
| 98 | do funit=10,100 |
|---|
| 99 | inquire(unit=funit, opened=is_used) |
|---|
| 100 | if (.not. is_used) exit |
|---|
| 101 | end do |
|---|
| 102 | open(funit,file='namelist.wps',status='old',form='formatted',err=1000) |
|---|
| 103 | read(funit,share) |
|---|
| 104 | read(funit,geogrid) |
|---|
| 105 | close(funit) |
|---|
| 106 | |
|---|
| 107 | dxkm = dx |
|---|
| 108 | dykm = dy |
|---|
| 109 | |
|---|
| 110 | known_lat = ref_lat |
|---|
| 111 | known_lon = ref_lon |
|---|
| 112 | known_x = ref_x |
|---|
| 113 | known_y = ref_y |
|---|
| 114 | |
|---|
| 115 | ! Convert wrf_core to uppercase letters |
|---|
| 116 | do i=1,3 |
|---|
| 117 | if (ichar(wrf_core(i:i)) >= 97) wrf_core(i:i) = char(ichar(wrf_core(i:i))-32) |
|---|
| 118 | end do |
|---|
| 119 | |
|---|
| 120 | ! Before doing anything else, we must have a valid grid type |
|---|
| 121 | gridtype = ' ' |
|---|
| 122 | if (wrf_core == 'ARW') then |
|---|
| 123 | gridtype = 'C' |
|---|
| 124 | dyn_opt = 2 |
|---|
| 125 | else if (wrf_core == 'NMM') then |
|---|
| 126 | gridtype = 'E' |
|---|
| 127 | dyn_opt = 4 |
|---|
| 128 | end if |
|---|
| 129 | |
|---|
| 130 | if (gridtype /= 'C' .and. gridtype /= 'E') then |
|---|
| 131 | write(6,*) 'A valid wrf_core must be specified in the namelist. '// & |
|---|
| 132 | 'Currently, only "ARW" and "NMM" are supported.' |
|---|
| 133 | stop |
|---|
| 134 | end if |
|---|
| 135 | |
|---|
| 136 | ! |
|---|
| 137 | ! Currently, plotgrids.exe does not work with NMM domains |
|---|
| 138 | ! |
|---|
| 139 | if (gridtype == 'E') then |
|---|
| 140 | write(6,*) |
|---|
| 141 | write(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' |
|---|
| 142 | write(6,*) '! This version of the plotgrids utility does not work for NMM domains. !' |
|---|
| 143 | write(6,*) '! An NMM-specific plotgrids may or may not be under development. !' |
|---|
| 144 | write(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' |
|---|
| 145 | write(6,*) |
|---|
| 146 | stop |
|---|
| 147 | end if |
|---|
| 148 | |
|---|
| 149 | if (max_dom > MAX_DOMAINS) then |
|---|
| 150 | write(6,*) 'In namelist, max_dom must be <= ',MAX_DOMAINS,'. To run with more'// & |
|---|
| 151 | ' than ',MAX_DOMAINS,' domains, increase the MAX_DOMAINS parameter.' |
|---|
| 152 | stop |
|---|
| 153 | end if |
|---|
| 154 | |
|---|
| 155 | ! Every domain must have a valid parent id |
|---|
| 156 | do i=2,max_dom |
|---|
| 157 | if (parent_id(i) <= 0 .or. parent_id(i) >= i) then |
|---|
| 158 | write(6,*) 'In namelist, the parent_id of domain ',i,' must be in '// & |
|---|
| 159 | 'the range 1 to ',i-1 |
|---|
| 160 | stop |
|---|
| 161 | end if |
|---|
| 162 | end do |
|---|
| 163 | |
|---|
| 164 | ! Convert map_proj to uppercase letters |
|---|
| 165 | do i=1,len(map_proj) |
|---|
| 166 | if (ichar(map_proj(i:i)) >= 97) map_proj(i:i) = char(ichar(map_proj(i:i))-32) |
|---|
| 167 | end do |
|---|
| 168 | |
|---|
| 169 | ! Assign parameters to module variables |
|---|
| 170 | if ((index(map_proj, 'LAMBERT') /= 0) .and. & |
|---|
| 171 | (len_trim(map_proj) == len('LAMBERT'))) then |
|---|
| 172 | iproj_type = PROJ_LC |
|---|
| 173 | rot=truelat1 |
|---|
| 174 | polat=truelat2 |
|---|
| 175 | jproj = 3 |
|---|
| 176 | |
|---|
| 177 | else if ((index(map_proj, 'MERCATOR') /= 0) .and. & |
|---|
| 178 | (len_trim(map_proj) == len('MERCATOR'))) then |
|---|
| 179 | iproj_type = PROJ_MERC |
|---|
| 180 | rot=0. |
|---|
| 181 | polat=0. |
|---|
| 182 | jproj = 9 |
|---|
| 183 | |
|---|
| 184 | else if ((index(map_proj, 'POLAR') /= 0) .and. & |
|---|
| 185 | (len_trim(map_proj) == len('POLAR'))) then |
|---|
| 186 | iproj_type = PROJ_PS |
|---|
| 187 | rot=0. |
|---|
| 188 | polat=SIGN(90., ref_lat) |
|---|
| 189 | jproj = 1 |
|---|
| 190 | |
|---|
| 191 | else if ((index(map_proj, 'ROTATED_LL') /= 0) .and. & |
|---|
| 192 | (len_trim(map_proj) == len('ROTATED_LL'))) then |
|---|
| 193 | iproj_type = PROJ_ROTLL |
|---|
| 194 | |
|---|
| 195 | else |
|---|
| 196 | write(6,*) 'In namelist, invalid map_proj specified. Valid '// & |
|---|
| 197 | 'projections are "lambert", "mercator", "polar", '// & |
|---|
| 198 | 'and "rotated_ll".' |
|---|
| 199 | end if |
|---|
| 200 | |
|---|
| 201 | n_domains = max_dom |
|---|
| 202 | |
|---|
| 203 | do i=1,n_domains |
|---|
| 204 | ixdim(i) = e_we(i) - s_we(i) + 1 |
|---|
| 205 | jydim(i) = e_sn(i) - s_sn(i) + 1 |
|---|
| 206 | end do |
|---|
| 207 | |
|---|
| 208 | if (gridtype == 'E') then |
|---|
| 209 | phi = dykm*real(jydim(1)-1)/2. |
|---|
| 210 | lambda = dxkm*real(ixdim(1)-1) |
|---|
| 211 | end if |
|---|
| 212 | |
|---|
| 213 | ! If the user hasn't supplied a known_x and known_y, assume the center of domain 1 |
|---|
| 214 | if (known_x == NAN) known_x = ixdim(1) / 2. |
|---|
| 215 | if (known_y == NAN) known_y = jydim(1) / 2. |
|---|
| 216 | |
|---|
| 217 | ! Checks specific to C grid |
|---|
| 218 | if (gridtype == 'C') then |
|---|
| 219 | |
|---|
| 220 | ! C grid does not support the rotated lat/lon projection |
|---|
| 221 | if (iproj_type == PROJ_ROTLL) then |
|---|
| 222 | write(6,*) 'Rotated lat/lon projection is not supported for the ARW core. '// & |
|---|
| 223 | 'Valid projecitons are "lambert", "mercator", and "polar".' |
|---|
| 224 | stop |
|---|
| 225 | end if |
|---|
| 226 | |
|---|
| 227 | ! Check that nests have an acceptable number of grid points in each dimension |
|---|
| 228 | do i=2,n_domains |
|---|
| 229 | rparent_gridpts = real(ixdim(i)-1)/real(parent_grid_ratio(i)) |
|---|
| 230 | if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then |
|---|
| 231 | write(6,*) 'For nest ',i,' (e_we-s_we+1) must be one greater than an '// & |
|---|
| 232 | 'interger multiple of the parent_grid_ratio.' |
|---|
| 233 | stop |
|---|
| 234 | end if |
|---|
| 235 | rparent_gridpts = real(jydim(i)-1)/real(parent_grid_ratio(i)) |
|---|
| 236 | if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then |
|---|
| 237 | write(6,*) 'For nest ',i,' (e_sn-s_sn+1) must be one greater than an '// & |
|---|
| 238 | 'interger multiple of the parent_grid_ratio.' |
|---|
| 239 | stop |
|---|
| 240 | end if |
|---|
| 241 | end do |
|---|
| 242 | end if |
|---|
| 243 | |
|---|
| 244 | do i=1,n_domains |
|---|
| 245 | parent_ll_x(i) = real(i_parent_start(i)) |
|---|
| 246 | parent_ll_y(i) = real(j_parent_start(i)) |
|---|
| 247 | parent_ur_x(i) = real(i_parent_start(i))+real(ixdim(i))/real(parent_grid_ratio(i))-1. |
|---|
| 248 | parent_ur_y(i) = real(j_parent_start(i))+real(jydim(i))/real(parent_grid_ratio(i))-1. |
|---|
| 249 | end do |
|---|
| 250 | |
|---|
| 251 | call map_init(map_projection) |
|---|
| 252 | |
|---|
| 253 | call map_set(iproj_type, map_projection, & |
|---|
| 254 | lat1=known_lat, & |
|---|
| 255 | lon1=known_lon, & |
|---|
| 256 | knowni=known_x, & |
|---|
| 257 | knownj=known_y, & |
|---|
| 258 | dx=dx, & |
|---|
| 259 | stdlon=stand_lon, & |
|---|
| 260 | truelat1=truelat1, & |
|---|
| 261 | truelat2=truelat2, & |
|---|
| 262 | ixdim=ixdim(1), & |
|---|
| 263 | jydim=jydim(1)) |
|---|
| 264 | |
|---|
| 265 | call ij_to_latlon(map_projection, 1., 1., rlats, rlons) |
|---|
| 266 | call ij_to_latlon(map_projection, real(e_we(1)), real(e_sn(1)) , rlate, rlone) |
|---|
| 267 | |
|---|
| 268 | call opngks |
|---|
| 269 | |
|---|
| 270 | ! Set some colors |
|---|
| 271 | call gscr(1, 0, 1.00, 1.00, 1.00) |
|---|
| 272 | call gscr(1, 1, 0.00, 0.00, 0.00) |
|---|
| 273 | |
|---|
| 274 | ! Do not grind them with details |
|---|
| 275 | jgrid=10 |
|---|
| 276 | jlts=-2 |
|---|
| 277 | iusout=1 |
|---|
| 278 | idot=0 |
|---|
| 279 | |
|---|
| 280 | call supmap(jproj,polat,stand_lon,rot,& |
|---|
| 281 | rlats,rlons,rlate,rlone, & |
|---|
| 282 | jlts,jgrid,iusout,idot,ier) |
|---|
| 283 | |
|---|
| 284 | call setusv('LW',1000) |
|---|
| 285 | call perim(e_we(1)-1,1,e_sn(1)-1,1) |
|---|
| 286 | call getset(xa,xb,ya,yb,xxa,xxy,yya,yyb,ltype) |
|---|
| 287 | call set (xa,xb,ya,yb, & |
|---|
| 288 | 1.,real(e_we(1)),1.,real(e_sn(1)),ltype) |
|---|
| 289 | |
|---|
| 290 | do idom = 2 , max_dom |
|---|
| 291 | call getxy ( xs, xe, ys, ye, & |
|---|
| 292 | idom , max_dom , & |
|---|
| 293 | e_we , e_sn , & |
|---|
| 294 | parent_id , parent_grid_ratio , & |
|---|
| 295 | i_parent_start , j_parent_start ) |
|---|
| 296 | call line ( xs , ys , xe , ys ) |
|---|
| 297 | call line ( xe , ys , xe , ye ) |
|---|
| 298 | call line ( xe , ye , xs , ye ) |
|---|
| 299 | call line ( xs , ye , xs , ys ) |
|---|
| 300 | end do |
|---|
| 301 | |
|---|
| 302 | call frame |
|---|
| 303 | |
|---|
| 304 | call clsgks |
|---|
| 305 | |
|---|
| 306 | stop |
|---|
| 307 | |
|---|
| 308 | 1000 write(6,*) 'Error opening namelist.wps' |
|---|
| 309 | stop |
|---|
| 310 | |
|---|
| 311 | end program plotgrids |
|---|
| 312 | |
|---|
| 313 | subroutine getxy ( xs, xe, ys, ye, & |
|---|
| 314 | dom_id , num_domains , & |
|---|
| 315 | e_we , e_sn , & |
|---|
| 316 | parent_id , parent_grid_ratio , & |
|---|
| 317 | i_parent_start , j_parent_start ) |
|---|
| 318 | |
|---|
| 319 | implicit none |
|---|
| 320 | |
|---|
| 321 | integer , intent(in) :: dom_id |
|---|
| 322 | integer , intent(in) :: num_domains |
|---|
| 323 | integer , intent(in) , dimension(num_domains):: e_we , e_sn , & |
|---|
| 324 | parent_id , parent_grid_ratio , & |
|---|
| 325 | i_parent_start , j_parent_start |
|---|
| 326 | real , intent(out) :: xs, xe, ys, ye |
|---|
| 327 | |
|---|
| 328 | |
|---|
| 329 | ! local vars |
|---|
| 330 | |
|---|
| 331 | integer :: idom |
|---|
| 332 | |
|---|
| 333 | xs = 0. |
|---|
| 334 | xe = e_we(dom_id) -1 |
|---|
| 335 | ys = 0. |
|---|
| 336 | ye = e_sn(dom_id) -1 |
|---|
| 337 | |
|---|
| 338 | idom = dom_id |
|---|
| 339 | compute_xy : DO |
|---|
| 340 | |
|---|
| 341 | xs = (i_parent_start(idom) + xs -1 ) / & |
|---|
| 342 | real(parent_grid_ratio(parent_id(idom))) |
|---|
| 343 | xe = xe / real(parent_grid_ratio(idom)) |
|---|
| 344 | |
|---|
| 345 | ys = (j_parent_start(idom) + ys -1 ) / & |
|---|
| 346 | real(parent_grid_ratio(parent_id(idom))) |
|---|
| 347 | ye = ye / real(parent_grid_ratio(idom)) |
|---|
| 348 | |
|---|
| 349 | idom = parent_id(idom) |
|---|
| 350 | if ( idom .EQ. 1 ) then |
|---|
| 351 | exit compute_xy |
|---|
| 352 | end if |
|---|
| 353 | |
|---|
| 354 | END DO compute_xy |
|---|
| 355 | |
|---|
| 356 | xs = xs + 1 |
|---|
| 357 | xe = xs + xe |
|---|
| 358 | ys = ys + 1 |
|---|
| 359 | ye = ys + ye |
|---|
| 360 | |
|---|
| 361 | end subroutine getxy |
|---|
