1 !***********************************************************************
2 !* GNU Lesser General Public License
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
6 !* FMS is free software: you can redistribute it and/or modify it under
7 !* the terms of the GNU Lesser General Public License as published by
8 !* the Free Software Foundation, either version 3 of the License, or (at
9 !* your option) any later version.
11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 !* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 !* You should have received a copy of the GNU Lesser General Public
17 !* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18 !***********************************************************************
19 !> @defgroup grid_mod grid_mod
21 !> @brief Routines for grid calculations
24 #ifdef use_deprecated_io
26 use mpp_mod, only : mpp_root_pe, uppercase, lowercase, FATAL, NOTE, mpp_error
27 use constants_mod, only : PI, radius
28 use fms_io_mod, only : get_great_circle_algorithm, get_global_att_value, string, &
29 field_exist, field_size, read_data
30 use mosaic_mod, only : get_mosaic_ntiles, get_mosaic_xgrid_size, get_mosaic_grid_sizes, &
31 get_mosaic_xgrid, calc_mosaic_grid_area, calc_mosaic_grid_great_circle_area
33 ! the following two use statement are only needed for define_cube_mosaic
34 use mpp_domains_mod, only : domain2d, mpp_define_mosaic, mpp_get_compute_domain, &
35 mpp_get_global_domain, domainUG, mpp_pass_SG_to_UG
36 use mosaic_mod, only : get_mosaic_ncontacts, get_mosaic_contact
40 ! ==== public interfaces =====================================================
41 ! grid dimension inquiry subroutines
42 public :: get_grid_ntiles
43 public :: get_grid_size
44 ! grid geometry inquiry subroutines
45 public :: get_grid_cell_centers
46 public :: get_grid_cell_vertices
47 ! grid area inquiry subroutines
48 public :: get_grid_cell_area
49 public :: get_grid_comp_area
50 ! decompose cubed sphere domains -- probably does not belong here, but it should
51 ! be in some place available for component models
52 public :: define_cube_mosaic
53 ! ==== end of public interfaces ==============================================
55 !> returns horizontal sizes of the grid
57 interface get_grid_size
58 module procedure get_grid_size_for_all_tiles
59 module procedure get_grid_size_for_one_tile
61 !> returns number of tiles
63 interface get_grid_cell_vertices
64 module procedure get_grid_cell_vertices_1D
65 module procedure get_grid_cell_vertices_2D
66 module procedure get_grid_cell_vertices_UG
70 interface get_grid_cell_centers
71 module procedure get_grid_cell_centers_1D
72 module procedure get_grid_cell_centers_2D
73 module procedure get_grid_cell_centers_UG
77 interface get_grid_cell_area
78 module procedure get_grid_cell_area_SG
79 module procedure get_grid_cell_area_UG
80 end interface get_grid_cell_area
83 interface get_grid_comp_area
84 module procedure get_grid_comp_area_SG
85 module procedure get_grid_comp_area_UG
86 end interface get_grid_comp_area
88 !> @addtogroup grid_mod
91 ! ==== module constants ======================================================
92 character(len=*), parameter :: &
93 module_name = 'grid_mod'
95 ! Include variable "version" to be written to log file.
96 #include<file_version.h>
98 character(len=*), parameter :: &
99 grid_dir = 'INPUT/', & !< root directory for all grid files
100 grid_file = 'INPUT/grid_spec.nc' !< name of the grid spec file
102 integer, parameter :: &
103 MAX_NAME = 256, & !< max length of the variable names
104 MAX_FILE = 1024, & !< max length of the file names
109 integer, parameter :: BUFSIZE = 1048576 !< This is used to control memory usage in get_grid_comp_area
110 !! We may change this to a namelist variable is needed.
112 ! ==== module variables ======================================================
113 integer :: grid_version = -1
114 logical :: great_circle_algorithm = .FALSE.
115 logical :: first_call = .TRUE.
120 function get_grid_version()
121 integer :: get_grid_version
124 great_circle_algorithm = get_great_circle_algorithm()
128 if(grid_version<0) then
129 if(field_exist(grid_file, 'geolon_t')) then
130 grid_version = VERSION_0
131 else if(field_exist(grid_file, 'x_T')) then
132 grid_version = VERSION_1
133 else if(field_exist(grid_file, 'ocn_mosaic_file') ) then
134 grid_version = VERSION_2
136 call mpp_error(FATAL, module_name//&
137 & '/get_grid_version: Can''t determine the version of the grid spec:'// &
138 & ' none of "x_T", "geolon_t", or "ocn_mosaic_file" exist in file "'//trim(grid_file)//'"')
141 get_grid_version = grid_version
142 end function get_grid_version
145 ! ============================================================================
146 ! ============================================================================
147 !> Returns number of tiles for a given component
148 subroutine get_grid_ntiles(component,ntiles)
149 character(len=*) :: component
150 integer, intent(out) :: ntiles
153 character(len=MAX_FILE) :: component_mosaic
155 select case (get_grid_version())
156 case(VERSION_0,VERSION_1)
159 call read_data(grid_file,trim(lowercase(component))//'_mosaic_file',component_mosaic)
160 ntiles = get_mosaic_ntiles(grid_dir//trim(component_mosaic))
162 end subroutine get_grid_ntiles
165 ! ============================================================================
166 ! ============================================================================
167 !> Returns size of the grid for each of the tiles
168 subroutine get_grid_size_for_all_tiles(component,nx,ny)
169 character(len=*) :: component
170 integer, intent(inout) :: nx(:),ny(:)
173 integer :: siz(4) ! for the size of external fields
174 character(len=MAX_NAME) :: varname1, varname2
175 character(len=MAX_FILE) :: component_mosaic
177 varname1 = 'AREA_'//trim(uppercase(component))
178 varname2 = trim(lowercase(component))//'_mosaic_file'
180 select case (get_grid_version())
181 case(VERSION_0,VERSION_1)
182 call field_size(grid_file, varname1, siz)
183 nx(1) = siz(1); ny(1)=siz(2)
184 case(VERSION_2) ! mosaic file
185 call read_data(grid_file,varname2, component_mosaic)
186 call get_mosaic_grid_sizes(grid_dir//trim(component_mosaic),nx,ny)
188 end subroutine get_grid_size_for_all_tiles
191 ! ============================================================================
192 ! ============================================================================
193 !> Returns size of the grid for one of the tiles
194 subroutine get_grid_size_for_one_tile(component,tile,nx,ny)
195 character(len=*) :: component
196 integer, intent(in) :: tile
197 integer, intent(inout) :: nx,ny
200 integer, allocatable :: nnx(:), nny(:)
203 call get_grid_ntiles(component, ntiles)
204 if(tile>0.and.tile<=ntiles) then
205 allocate(nnx(ntiles),nny(ntiles))
206 call get_grid_size_for_all_tiles(component,nnx,nny)
207 nx = nnx(tile); ny = nny(tile)
210 call mpp_error(FATAL, 'get_grid_size: requested tile index '// &
211 & trim(string(tile))//' is out of bounds (1:'//trim(string(ntiles))//')')
213 end subroutine get_grid_size_for_one_tile
215 ! ============================================================================
216 ! ============================================================================
217 !> Return grid cell area for the specified model component and tile
218 subroutine get_grid_cell_area_SG(component, tile, cellarea, domain)
219 character(len=*), intent(in) :: component
220 integer , intent(in) :: tile
221 real , intent(inout) :: cellarea(:,:)
222 type(domain2d) , intent(in), optional :: domain
225 integer :: nlon, nlat
226 real, allocatable :: glonb(:,:), glatb(:,:)
228 select case(get_grid_version())
229 case(VERSION_0,VERSION_1)
230 select case(trim(component))
232 call read_data(grid_file, 'AREA_LND_CELL', cellarea, &
233 no_domain=.not.present(domain), domain=domain)
235 call read_data(grid_file, 'AREA_'//trim(uppercase(component)),cellarea,&
236 no_domain=.not.present(domain),domain=domain)
238 call mpp_error(FATAL, module_name//'/get_grid_cell_area: Illegal component name "'//trim(component) &
239 & //'": must be one of ATM, LND, or OCN')
242 cellarea = cellarea*4.*PI*radius**2
244 if (present(domain)) then
245 call mpp_get_compute_domain(domain,xsize=nlon,ysize=nlat)
247 call get_grid_size(component,tile,nlon,nlat)
249 allocate(glonb(nlon+1,nlat+1),glatb(nlon+1,nlat+1))
250 call get_grid_cell_vertices(component, tile, glonb, glatb, domain)
251 if (great_circle_algorithm) then
252 call calc_mosaic_grid_great_circle_area(glonb*pi/180.0, glatb*pi/180.0, cellarea)
254 call calc_mosaic_grid_area(glonb*pi/180.0, glatb*pi/180.0, cellarea)
256 deallocate(glonb,glatb)
259 end subroutine get_grid_cell_area_SG
261 ! ============================================================================
262 ! ============================================================================
263 !> Get the area of the component per grid cell
264 subroutine get_grid_comp_area_SG(component,tile,area,domain)
265 character(len=*) :: component
266 integer, intent(in) :: tile
267 real, intent(inout) :: area(:,:)
268 type(domain2d), intent(in), optional :: domain
270 integer :: n_xgrid_files ! number of exchange grid files in the mosaic
271 integer :: siz(4), nxgrid
273 integer, allocatable :: i1(:), j1(:), i2(:), j2(:)
274 real, allocatable :: xgrid_area(:)
275 real, allocatable :: rmask(:,:)
276 character(len=MAX_NAME) :: &
277 xgrid_name, & ! name of the variable holding xgrid names
278 tile_name, & ! name of the tile
279 xgrid_file, & ! name of the current xgrid file
280 mosaic_name,& ! name of the mosaic
283 character(len=4096) :: attvalue
284 character(len=MAX_NAME), allocatable :: nest_tile_name(:)
285 integer :: is,ie,js,je ! boundaries of our domain
286 integer :: i0, j0 ! offsets for x and y, respectively
287 integer :: num_nest_tile, ntiles
289 integer :: found_xgrid_files ! how many xgrid files we actually found in the grid spec
290 integer :: ibegin, iend, bsize, l
292 select case (get_grid_version())
293 case(VERSION_0,VERSION_1)
294 select case(component)
296 call read_data(grid_file,'AREA_ATM',area, no_domain=.not.present(domain),domain=domain)
298 allocate(rmask(size(area,1),size(area,2)))
299 call read_data(grid_file,'AREA_OCN',area, no_domain=.not.present(domain),domain=domain)
300 call read_data(grid_file,'wet', rmask,no_domain=.not.present(domain),domain=domain)
304 call read_data(grid_file,'AREA_LND',area,no_domain=.not.present(domain),domain=domain)
306 call mpp_error(FATAL, module_name// &
307 & '/get_grid_comp_area: Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN')
309 case(VERSION_2) ! mosaic gridspec
310 select case (component)
312 ! just read the grid cell area and return
313 call get_grid_cell_area(component,tile,area)
316 xgrid_name = 'aXl_file'
317 call read_data(grid_file, 'lnd_mosaic', mosaic_name)
318 tile_name = trim(mosaic_name)//'_tile'//char(tile+ichar('0'))
320 xgrid_name = 'aXo_file'
321 call read_data(grid_file, 'ocn_mosaic', mosaic_name)
322 tile_name = trim(mosaic_name)//'_tile'//char(tile+ichar('0'))
324 call mpp_error(FATAL, module_name// &
325 & '/get_grid_comp_area: Illegal component name "'//trim(component)//'": must be one of ATM, LND, or OCN')
327 ! get the boundaries of the requested domain
328 if(present(domain)) then
329 call mpp_get_compute_domain(domain,is,ie,js,je)
332 call get_grid_size(component,tile,ie,je)
336 if (size(area,1)/=ie-is+1.or.size(area,2)/=je-js+1) &
337 call mpp_error(FATAL, module_name// &
338 & '/get_grid_comp_area: size of the output argument "area" is not consistent with the domain')
341 call read_data(grid_file, 'atm_mosaic', mosaic_name)
342 call read_data(grid_file,'atm_mosaic_file',mosaic_file)
343 mosaic_file = grid_dir//trim(mosaic_file)
344 ntiles = get_mosaic_ntiles(trim(mosaic_file))
345 allocate(nest_tile_name(ntiles))
348 call read_data(mosaic_file, 'gridfiles', tilefile, level=n)
349 tilefile = grid_dir//trim(tilefile)
350 if( get_global_att_value(tilefile, "nest_grid", attvalue) ) then
351 if(trim(attvalue) == "TRUE") then
352 num_nest_tile = num_nest_tile + 1
353 nest_tile_name(num_nest_tile) = trim(mosaic_name)//'_tile'//char(n+ichar('0'))
354 else if(trim(attvalue) .NE. "FALSE") then
355 call mpp_error(FATAL, module_name//'/get_grid_comp_area: value of global attribute nest_grid in file'// &
356 & trim(tilefile)//' should be TRUE of FALSE')
361 if(field_exist(grid_file,xgrid_name)) then
362 ! get the number of the exchange-grid files
363 call field_size(grid_file,xgrid_name,siz)
364 n_xgrid_files = siz(2)
365 found_xgrid_files = 0
366 ! loop through all exchange grid files
367 do n = 1, n_xgrid_files
368 ! get the name of the current exchange grid file
369 call read_data(grid_file,xgrid_name,xgrid_file,level=n)
370 ! skip the rest of the loop if the name of the current tile isn't found
371 ! in the file name, but check this only if there is more than 1 tile
372 if(n_xgrid_files>1) then
373 if(index(xgrid_file,trim(tile_name))==0) cycle
375 found_xgrid_files = found_xgrid_files + 1
376 !---make sure the atmosphere grid is not a nested grid
378 do m = 1, num_nest_tile
379 if(index(xgrid_file, trim(nest_tile_name(m))) .NE. 0) then
386 ! finally read the exchange grid
387 nxgrid = get_mosaic_xgrid_size(grid_dir//xgrid_file)
388 if(nxgrid < BUFSIZE) then
389 allocate(i1(nxgrid), j1(nxgrid), i2(nxgrid), j2(nxgrid), xgrid_area(nxgrid))
391 allocate(i1(BUFSIZE), j1(BUFSIZE), i2(BUFSIZE), j2(BUFSIZE), xgrid_area(BUFSIZE))
394 do l = 1,nxgrid,BUFSIZE
395 bsize = min(BUFSIZE, nxgrid-l+1)
396 iend = ibegin + bsize - 1
397 call get_mosaic_xgrid(grid_dir//xgrid_file, i1(1:bsize), j1(1:bsize), i2(1:bsize), j2(1:bsize), &
398 xgrid_area(1:bsize), ibegin, iend)
399 ! and sum the exchange grid areas
402 if (i<is.or.i>ie) cycle
403 if (j<js.or.j>je) cycle
404 area(i+i0,j+j0) = area(i+i0,j+j0) + xgrid_area(m)
408 deallocate(i1, j1, i2, j2, xgrid_area)
410 if (found_xgrid_files == 0) &
411 call mpp_error(FATAL, 'get_grid_comp_area: no xgrid files were found for component '// &
412 & trim(component)//' (mosaic name is '//trim(mosaic_name)//')')
415 deallocate(nest_tile_name)
418 area = area*4.*PI*radius**2
419 end subroutine get_grid_comp_area_SG
421 !======================================================================
422 subroutine get_grid_cell_area_UG(component, tile, cellarea, SG_domain, UG_domain)
423 character(len=*), intent(in) :: component
424 integer , intent(in) :: tile
425 real , intent(inout) :: cellarea(:)
426 type(domain2d) , intent(in) :: SG_domain
427 type(domainUG) , intent(in) :: UG_domain
428 integer :: is, ie, js, je
429 real, allocatable :: SG_area(:,:)
431 call mpp_get_compute_domain(SG_domain, is, ie, js, je)
432 allocate(SG_area(is:ie, js:je))
433 call get_grid_cell_area_SG(component, tile, SG_area, SG_domain)
434 call mpp_pass_SG_to_UG(UG_domain, SG_area, cellarea)
437 end subroutine get_grid_cell_area_UG
439 subroutine get_grid_comp_area_UG(component, tile, area, SG_domain, UG_domain)
440 character(len=*), intent(in) :: component
441 integer , intent(in) :: tile
442 real , intent(inout) :: area(:)
443 type(domain2d) , intent(in) :: SG_domain
444 type(domainUG) , intent(in) :: UG_domain
445 integer :: is, ie, js, je
446 real, allocatable :: SG_area(:,:)
448 call mpp_get_compute_domain(SG_domain, is, ie, js, je)
449 allocate(SG_area(is:ie, js:je))
450 call get_grid_comp_area_SG(component, tile, SG_area, SG_domain)
451 call mpp_pass_SG_to_UG(UG_domain, SG_area, area)
454 end subroutine get_grid_comp_area_UG
457 ! ============================================================================
458 ! ============================================================================
459 !> Returns arrays of global grid cell boundaries for given model component and
460 !! mosaic tile number.
462 !> @note In the case of non-lat-lon grid the returned coordinates may have be not so
463 !! meaningful, by the very nature of such grids. But presumably these 1D coordinate
464 !! arrays are good enough for diag axis and such.
465 subroutine get_grid_cell_vertices_1D(component, tile, glonb, glatb)
466 character(len=*), intent(in) :: component
467 integer, intent(in) :: tile
468 real, intent(inout) :: glonb(:),glatb(:)
470 integer :: nlon, nlat
471 integer :: start(4), nread(4)
472 real, allocatable :: tmp(:,:), x_vert_t(:,:,:), y_vert_t(:,:,:)
473 character(len=MAX_FILE) :: filename1, filename2
475 call get_grid_size_for_one_tile(component, tile, nlon, nlat)
476 if (size(glonb(:))/=nlon+1) &
477 call mpp_error (FATAL, module_name// &
478 & '/get_grid_cell_vertices_1D: Size of argument "glonb" is not consistent with the grid size')
479 if (size(glatb(:))/=nlat+1) &
480 call mpp_error (FATAL, module_name// &
481 & '/get_grid_cell_vertices_1D: Size of argument "glatb" is not consistent with the grid size')
482 if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
483 call mpp_error(FATAL, module_name//'/get_grid_cell_vertices_1D: Illegal component name "'// &
484 & trim(component)//'": must be one of ATM, LND, or OCN')
487 select case(get_grid_version())
489 select case(trim(component))
491 call read_data(grid_file, 'xb'//lowercase(component(1:1)), glonb, no_domain=.true.)
492 call read_data(grid_file, 'yb'//lowercase(component(1:1)), glatb, no_domain=.true.)
494 call read_data(grid_file, "gridlon_vert_t", glonb, no_domain=.true.)
495 call read_data(grid_file, "gridlat_vert_t", glatb, no_domain=.true.)
498 select case(trim(component))
500 call read_data(grid_file, 'xb'//lowercase(component(1:1)), glonb, no_domain=.true.)
501 call read_data(grid_file, 'yb'//lowercase(component(1:1)), glatb, no_domain=.true.)
503 allocate (x_vert_t(nlon,1,2), y_vert_t(1,nlat,2) )
505 nread(1) = nlon; nread(2) = 1; start(3) = 1
506 call read_data(grid_file, "x_vert_T", x_vert_t(:,:,1), start, nread, no_domain=.TRUE.)
507 nread(1) = nlon; nread(2) = 1; start(3) = 2
508 call read_data(grid_file, "x_vert_T", x_vert_t(:,:,2), start, nread, no_domain=.TRUE.)
510 nread(1) = 1; nread(2) = nlat; start(3) = 1
511 call read_data(grid_file, "y_vert_T", y_vert_t(:,:,1), start, nread, no_domain=.TRUE.)
512 nread(1) = 1; nread(2) = nlat; start(3) = 4
513 call read_data(grid_file, "y_vert_T", y_vert_t(:,:,2), start, nread, no_domain=.TRUE.)
514 glonb(1:nlon) = x_vert_t(1:nlon,1,1)
515 glonb(nlon+1) = x_vert_t(nlon,1,2)
516 glatb(1:nlat) = y_vert_t(1,1:nlat,1)
517 glatb(nlat+1) = y_vert_t(1,nlat,2)
518 deallocate(x_vert_t, y_vert_t)
521 ! get the name of the mosaic file for the component
522 call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
523 filename1=grid_dir//trim(filename1)
524 ! get the name of the grid file for the component and tile
525 call read_data(filename1, 'gridfiles', filename2, level=tile)
526 filename2 = grid_dir//trim(filename2)
530 allocate( tmp(2*nlon+1,1) )
531 call read_data(filename2, "x", tmp, start, nread, no_domain=.TRUE.)
532 glonb(1:nlon+1) = tmp(1:2*nlon+1:2,1)
534 allocate(tmp(1,2*nlat+1))
538 call read_data(filename2, "y", tmp, start, nread, no_domain=.TRUE.)
539 glatb(1:nlat+1) = tmp(1,1:2*nlat+1:2)
543 end subroutine get_grid_cell_vertices_1D
545 ! ============================================================================
546 ! ============================================================================
547 !> Returns cell vertices for the specified model component and mosaic tile number
548 subroutine get_grid_cell_vertices_2D(component, tile, lonb, latb, domain)
549 character(len=*), intent(in) :: component
550 integer, intent(in) :: tile
551 real, intent(inout) :: lonb(:,:),latb(:,:)
552 type(domain2d), optional, intent(in) :: domain
555 character(len=MAX_FILE) :: filename1, filename2
556 integer :: nlon, nlat
558 real, allocatable :: buffer(:), tmp(:,:), x_vert_t(:,:,:), y_vert_t(:,:,:)
559 integer :: is,ie,js,je ! boundaries of our domain
560 integer :: i0,j0 ! offsets for coordinates
562 integer :: start(4), nread(4)
564 call get_grid_size_for_one_tile(component, tile, nlon, nlat)
565 if (present(domain)) then
566 call mpp_get_compute_domain(domain,is,ie,js,je)
570 !--- domain normally should be present
571 call mpp_error (NOTE, module_name//'/get_grid_cell_vertices: domain is not present, global data will be read')
573 i0 = -is+1; j0 = -js+1
575 ! verify that lonb and latb sizes are consistent with the size of domain
576 if (size(lonb,1)/=ie-is+2.or.size(lonb,2)/=je-js+2) &
577 call mpp_error (FATAL, module_name// &
578 & '/get_grid_cell_vertices: Size of argument "lonb" is not consistent with the domain size')
579 if (size(latb,1)/=ie-is+2.or.size(latb,2)/=je-js+2) &
580 call mpp_error (FATAL, module_name// &
581 & '/get_grid_cell_vertices: Size of argument "latb" is not consistent with the domain size')
582 if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
583 call mpp_error(FATAL, module_name//'/get_grid_cell_vertices: Illegal component name "'// &
584 & trim(component)//'": must be one of ATM, LND, or OCN')
587 select case(get_grid_version())
589 select case(component)
591 allocate(buffer(max(nlon,nlat)+1))
592 ! read coordinates of grid cell vertices
593 call read_data(grid_file, 'xb'//lowercase(component(1:1)), buffer(1:nlon+1), no_domain=.true.)
596 lonb(i+i0,j+j0) = buffer(i)
599 call read_data(grid_file, 'yb'//lowercase(component(1:1)), buffer(1:nlat+1), no_domain=.true.)
602 latb(i+i0,j+j0) = buffer(j)
607 if (present(domain)) then
609 start(1) = is; start(2) = js
610 nread(1) = ie-is+2; nread(2) = je-js+2
611 call read_data(grid_file, 'geolon_vert_t', lonb, start, nread, no_domain=.true. )
612 call read_data(grid_file, 'geolat_vert_t', latb, start, nread, no_domain=.true. )
614 call read_data(grid_file, 'geolon_vert_t', lonb, no_domain=.TRUE. )
615 call read_data(grid_file, 'geolat_vert_t', latb, no_domain=.TRUE. )
619 select case(component)
621 allocate(buffer(max(nlon,nlat)+1))
622 ! read coordinates of grid cell vertices
623 call read_data(grid_file, 'xb'//lowercase(component(1:1)), buffer(1:nlon+1), no_domain=.true.)
626 lonb(i+i0,j+j0) = buffer(i)
629 call read_data(grid_file, 'yb'//lowercase(component(1:1)), buffer(1:nlat+1), no_domain=.true.)
632 latb(i+i0,j+j0) = buffer(j)
637 nlon=ie-is+1; nlat=je-js+1
638 allocate (x_vert_t(nlon,nlat,4), y_vert_t(nlon,nlat,4) )
639 call read_data(grid_file, 'x_vert_T', x_vert_t, no_domain=.not.present(domain), domain=domain )
640 call read_data(grid_file, 'y_vert_T', y_vert_t, no_domain=.not.present(domain), domain=domain )
641 lonb(1:nlon,1:nlat) = x_vert_t(1:nlon,1:nlat,1)
642 lonb(nlon+1,1:nlat) = x_vert_t(nlon,1:nlat,2)
643 lonb(1:nlon,nlat+1) = x_vert_t(1:nlon,nlat,4)
644 lonb(nlon+1,nlat+1) = x_vert_t(nlon,nlat,3)
645 latb(1:nlon,1:nlat) = y_vert_t(1:nlon,1:nlat,1)
646 latb(nlon+1,1:nlat) = y_vert_t(nlon,1:nlat,2)
647 latb(1:nlon,nlat+1) = y_vert_t(1:nlon,nlat,4)
648 latb(nlon+1,nlat+1) = y_vert_t(nlon,nlat,3)
649 deallocate(x_vert_t, y_vert_t)
652 ! get the name of the mosaic file for the component
653 call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
654 filename1=grid_dir//trim(filename1)
655 ! get the name of the grid file for the component and tile
656 call read_data(filename1, 'gridfiles', filename2, level=tile)
657 filename2 = grid_dir//trim(filename2)
658 if(PRESENT(domain)) then
659 call mpp_get_global_domain(domain, xbegin=isg, ybegin=jsg)
661 start(1) = 2*(is-isg+1) - 1; nread(1) = 2*(ie-is)+3
662 start(2) = 2*(js-jsg+1) - 1; nread(2) = 2*(je-js)+3
663 allocate(tmp(nread(1), nread(2)) )
664 call read_data(filename2, 'x', tmp, start, nread, no_domain=.TRUE.)
667 lonb(i,j) = tmp(2*i-1,2*j-1)
670 call read_data(filename2, 'y', tmp, start, nread, no_domain=.TRUE.)
673 latb(i,j) = tmp(2*i-1,2*j-1)
677 allocate(tmp(2*nlon+1,2*nlat+1))
678 call read_data(filename2, 'x', tmp, no_domain=.TRUE.)
681 lonb(i+i0,j+j0) = tmp(2*i-1,2*j-1)
684 call read_data(filename2, 'y', tmp, no_domain=.TRUE.)
687 latb(i+i0,j+j0) = tmp(2*i-1,2*j-1)
694 end subroutine get_grid_cell_vertices_2D
697 subroutine get_grid_cell_vertices_UG(component, tile, lonb, latb, SG_domain, UG_domain)
698 character(len=*), intent(in) :: component
699 integer, intent(in) :: tile
700 real, intent(inout) :: lonb(:,:),latb(:,:) ! The second dimension is 4
701 type(domain2d) , intent(in) :: SG_domain
702 type(domainUG) , intent(in) :: UG_domain
703 integer :: is, ie, js, je, i, j
704 real, allocatable :: SG_lonb(:,:), SG_latb(:,:), tmp(:,:,:)
706 call mpp_get_compute_domain(SG_domain, is, ie, js, je)
707 allocate(SG_lonb(is:ie+1, js:je+1))
708 allocate(SG_latb(is:ie+1, js:je+1))
709 allocate(tmp(is:ie,js:je,4))
710 call get_grid_cell_vertices_2D(component, tile, SG_lonb, SG_latb, SG_domain)
713 tmp(i,j,1) = SG_lonb(i,j)
714 tmp(i,j,2) = SG_lonb(i+1,j)
715 tmp(i,j,3) = SG_lonb(i+1,j+1)
716 tmp(i,j,4) = SG_lonb(i,j+1)
719 call mpp_pass_SG_to_UG(UG_domain, tmp, lonb)
722 tmp(i,j,1) = SG_latb(i,j)
723 tmp(i,j,2) = SG_latb(i+1,j)
724 tmp(i,j,3) = SG_latb(i+1,j+1)
725 tmp(i,j,4) = SG_latb(i,j+1)
728 call mpp_pass_SG_to_UG(UG_domain, tmp, latb)
731 deallocate(SG_lonb, SG_latb, tmp)
733 end subroutine get_grid_cell_vertices_UG
735 ! ============================================================================
736 !> Returns global coordinate arrays fro given model component and mosaic tile number
737 !! @note In the case of non-lat-lon grid those coordinates may have be not so
738 !! meaningful, by the very nature of such grids. But presumably these 1D coordinate
739 !! arrays are good enough for diag axis and such.
740 subroutine get_grid_cell_centers_1D(component, tile, glon, glat)
741 character(len=*), intent(in) :: component
742 integer, intent(in) :: tile
743 real, intent(inout) :: glon(:),glat(:)
744 integer :: nlon, nlat
745 integer :: start(4), nread(4)
746 real, allocatable :: tmp(:,:)
747 character(len=MAX_FILE) :: filename1, filename2
749 call get_grid_size_for_one_tile(component, tile, nlon, nlat)
750 if (size(glon(:))/=nlon) &
751 call mpp_error (FATAL, module_name// &
752 & '/get_grid_cell_centers_1D: Size of argument "glon" is not consistent with the grid size')
753 if (size(glat(:))/=nlat) &
754 call mpp_error (FATAL, module_name// &
755 & '/get_grid_cell_centers_1D: Size of argument "glat" is not consistent with the grid size')
756 if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
757 call mpp_error(FATAL, module_name//'/get_grid_cell_centers_1D: Illegal component name "'// &
758 & trim(component)//'": must be one of ATM, LND, or OCN')
761 select case(get_grid_version())
763 select case(trim(component))
765 call read_data(grid_file, 'xt'//lowercase(component(1:1)), glon, no_domain=.true.)
766 call read_data(grid_file, 'yt'//lowercase(component(1:1)), glat, no_domain=.true.)
768 call read_data(grid_file, "gridlon_t", glon, no_domain=.true.)
769 call read_data(grid_file, "gridlat_t", glat, no_domain=.true.)
772 select case(trim(component))
774 call read_data(grid_file, 'xt'//lowercase(component(1:1)), glon, no_domain=.true.)
775 call read_data(grid_file, 'yt'//lowercase(component(1:1)), glat, no_domain=.true.)
777 call read_data(grid_file, "grid_x_T", glon, no_domain=.true.)
778 call read_data(grid_file, "grid_y_T", glat, no_domain=.true.)
781 ! get the name of the mosaic file for the component
782 call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
783 filename1=grid_dir//trim(filename1)
784 ! get the name of the grid file for the component and tile
785 call read_data(filename1, 'gridfiles', filename2, level=tile)
786 filename2 = grid_dir//trim(filename2)
789 nread(1) = 2*nlon+1; start(2) = 2
790 allocate( tmp(2*nlon+1,1) )
791 call read_data(filename2, "x", tmp, start, nread, no_domain=.TRUE.)
792 glon(1:nlon) = tmp(2:2*nlon:2,1)
794 allocate(tmp(1, 2*nlat+1))
797 nread(2) = 2*nlat+1; start(1) = 2
798 call read_data(filename2, "y", tmp, start, nread, no_domain=.TRUE.)
799 glat(1:nlat) = tmp(1,2:2*nlat:2)
804 end subroutine get_grid_cell_centers_1D
806 ! ============================================================================
807 ! ============================================================================
808 !> Returns grid cell centers for specified model component and mosaic tile number
809 subroutine get_grid_cell_centers_2D(component, tile, lon, lat, domain)
810 character(len=*), intent(in) :: component
811 integer, intent(in) :: tile
812 real, intent(inout) :: lon(:,:),lat(:,:)
813 type(domain2d), intent(in), optional :: domain
815 character(len=MAX_FILE) :: filename1, filename2
816 integer :: nlon, nlat
818 real, allocatable :: buffer(:),tmp(:,:)
819 integer :: is,ie,js,je ! boundaries of our domain
820 integer :: i0,j0 ! offsets for coordinates
822 integer :: start(4), nread(4)
824 call get_grid_size_for_one_tile(component, tile, nlon, nlat)
825 if (present(domain)) then
826 call mpp_get_compute_domain(domain,is,ie,js,je)
830 !--- domain normally should be present
831 call mpp_error (NOTE, module_name//'/get_grid_cell_centers: domain is not present, global data will be read')
833 i0 = -is+1; j0 = -js+1
835 ! verify that lon and lat sizes are consistent with the size of domain
836 if (size(lon,1)/=ie-is+1.or.size(lon,2)/=je-js+1) &
837 call mpp_error (FATAL, module_name// &
838 & '/get_grid_cell_centers: Size of array "lon" is not consistent with the domain size')
839 if (size(lat,1)/=ie-is+1.or.size(lat,2)/=je-js+1) &
840 call mpp_error (FATAL, module_name// &
841 & '/get_grid_cell_centers: Size of array "lat" is not consistent with the domain size')
842 if(trim(component) .NE. 'ATM' .AND. component .NE. 'LND' .AND. component .NE. 'OCN') then
843 call mpp_error(FATAL, module_name//'/get_grid_cell_vertices: Illegal component name "'// &
844 & trim(component)//'": must be one of ATM, LND, or OCN')
847 select case(get_grid_version())
849 select case (trim(component))
851 allocate(buffer(max(nlon,nlat)))
852 ! read coordinates of grid cell vertices
853 call read_data(grid_file, 'xt'//lowercase(component(1:1)), buffer(1:nlon), no_domain=.true.)
856 lon(i+i0,j+j0) = buffer(i)
859 call read_data(grid_file, 'yt'//lowercase(component(1:1)), buffer(1:nlat), no_domain=.true.)
862 lat(i+i0,j+j0) = buffer(j)
867 call read_data(grid_file, 'geolon_t', lon, no_domain=.not.present(domain), domain=domain )
868 call read_data(grid_file, 'geolat_t', lat, no_domain=.not.present(domain), domain=domain )
871 select case(trim(component))
873 allocate(buffer(max(nlon,nlat)))
874 ! read coordinates of grid cell vertices
875 call read_data(grid_file, 'xt'//lowercase(component(1:1)), buffer(1:nlon), no_domain=.true.)
878 lon(i+i0,j+j0) = buffer(i)
881 call read_data(grid_file, 'yt'//lowercase(component(1:1)), buffer(1:nlat), no_domain=.true.)
884 lat(i+i0,j+j0) = buffer(j)
889 call read_data(grid_file, 'x_T', lon, no_domain=.not.present(domain), domain=domain )
890 call read_data(grid_file, 'y_T', lat, no_domain=.not.present(domain), domain=domain )
892 case(VERSION_2) ! mosaic grid file
893 ! get the name of the mosaic file for the component
894 call read_data(grid_file, trim(lowercase(component))//'_mosaic_file', filename1)
895 filename1=grid_dir//trim(filename1)
896 ! get the name of the grid file for the component and tile
897 call read_data(filename1, 'gridfiles', filename2, level=tile)
898 filename2 = grid_dir//trim(filename2)
899 if(PRESENT(domain)) then
900 call mpp_get_global_domain(domain, xbegin=isg, ybegin=jsg)
902 start(1) = 2*(is-isg+1) - 1; nread(1) = 2*(ie-is)+3
903 start(2) = 2*(js-jsg+1) - 1; nread(2) = 2*(je-js)+3
904 allocate(tmp(nread(1), nread(2)))
905 call read_data(filename2, 'x', tmp, start, nread, no_domain=.TRUE.)
908 lon(i,j) = tmp(2*i,2*j)
911 call read_data(filename2, 'y', tmp, start, nread, no_domain=.TRUE.)
914 lat(i,j) = tmp(2*i,2*j)
918 allocate(tmp(2*nlon+1,2*nlat+1))
919 call read_data(filename2, 'x', tmp, no_domain=.TRUE.)
922 lon(i+i0,j+j0) = tmp(2*i,2*j)
925 call read_data(filename2, 'y', tmp, no_domain=.TRUE.)
928 lat(i+i0,j+j0) = tmp(2*i,2*j)
935 end subroutine get_grid_cell_centers_2D
937 subroutine get_grid_cell_centers_UG(component, tile, lon, lat, SG_domain, UG_domain)
938 character(len=*), intent(in) :: component
939 integer, intent(in) :: tile
940 real, intent(inout) :: lon(:),lat(:)
941 type(domain2d) , intent(in) :: SG_domain
942 type(domainUG) , intent(in) :: UG_domain
943 integer :: is, ie, js, je
944 real, allocatable :: SG_lon(:,:), SG_lat(:,:)
946 call mpp_get_compute_domain(SG_domain, is, ie, js, je)
947 allocate(SG_lon(is:ie, js:je))
948 allocate(SG_lat(is:ie, js:je))
949 call get_grid_cell_centers_2D(component, tile, SG_lon, SG_lat, SG_domain)
950 call mpp_pass_SG_to_UG(UG_domain, SG_lon, lon)
951 call mpp_pass_SG_to_UG(UG_domain, SG_lat, lat)
952 deallocate(SG_lon, SG_lat)
954 end subroutine get_grid_cell_centers_UG
956 ! ============================================================================
957 ! ============================================================================
958 ! this subroutine probably does not belong in the grid_mod
959 !> Given a model component, a layout, and (optionally) a halo size, returns a
960 !! domain for current processor
961 subroutine define_cube_mosaic ( component, domain, layout, halo, maskmap )
962 character(len=*) , intent(in) :: component
963 type(domain2d) , intent(inout) :: domain
964 integer , intent(in) :: layout(2)
965 integer, optional, intent(in) :: halo
966 logical, optional, intent(in) :: maskmap(:,:,:)
968 ! ---- local constants
971 character(len=MAX_NAME) :: varname
972 character(len=MAX_FILE + len(grid_dir)) :: mosaic_file
973 integer :: ntiles ! number of tiles
974 integer :: ncontacts ! number of contacts between mosaic tiles
976 integer :: ng, pe_pos, npes ! halo size
977 integer, allocatable :: nlon(:), nlat(:), global_indices(:,:)
978 integer, allocatable :: pe_start(:), pe_end(:), layout_2d(:,:)
979 integer, allocatable :: tile1(:),tile2(:)
980 integer, allocatable :: is1(:),ie1(:),js1(:),je1(:)
981 integer, allocatable :: is2(:),ie2(:),js2(:),je2(:)
983 call get_grid_ntiles(component,ntiles)
984 allocate(nlon(ntiles), nlat(ntiles))
985 allocate(global_indices(4,ntiles))
986 allocate(pe_start(ntiles),pe_end(ntiles))
987 allocate(layout_2d(2,ntiles))
988 call get_grid_size(component,nlon,nlat)
990 pe_pos = mpp_root_pe()
992 global_indices(:,n) = (/ 1, nlon(n), 1, nlat(n) /)
993 layout_2d (:,n) = layout
994 if(present(maskmap)) then
995 npes = count(maskmap(:,:,n))
997 npes = layout(1)*layout(2)
1000 pe_end (n) = pe_pos + npes - 1
1001 pe_pos = pe_end(n) + 1
1004 varname=trim(lowercase(component))//'_mosaic_file'
1005 call read_data(grid_file,varname,mosaic_file(1:MAX_FILE))
1006 mosaic_file = grid_dir//mosaic_file(1:MAX_FILE)
1008 ! get the contact information from mosaic file
1009 ncontacts = get_mosaic_ncontacts(mosaic_file)
1010 allocate(tile1(ncontacts),tile2(ncontacts))
1011 allocate(is1(ncontacts),ie1(ncontacts),js1(ncontacts),je1(ncontacts))
1012 allocate(is2(ncontacts),ie2(ncontacts),js2(ncontacts),je2(ncontacts))
1013 call get_mosaic_contact(mosaic_file, tile1, tile2, &
1014 is1, ie1, js1, je1, is2, ie2, js2, je2)
1017 if(present(halo)) ng = halo
1018 ! create the domain2d variable
1019 call mpp_define_mosaic ( global_indices, layout_2d, domain, &
1020 ntiles, ncontacts, tile1, tile2, &
1021 is1, ie1, js1, je1, &
1022 is2, ie2, js2, je2, &
1023 pe_start=pe_start, pe_end=pe_end, symmetry=.true., &
1024 shalo = ng, nhalo = ng, whalo = ng, ehalo = ng, &
1025 maskmap = maskmap, &
1026 name = trim(component)//'Cubic-Sphere Grid' )
1028 deallocate(nlon,nlat,global_indices,pe_start,pe_end,layout_2d)
1029 deallocate(tile1,tile2)
1030 deallocate(is1,ie1,js1,je1)
1031 deallocate(is2,ie2,js2,je2)
1033 end subroutine define_cube_mosaic
1037 ! close documentation grouping