regrid_sfc  1.13.0
grids_IO.F90
Go to the documentation of this file.
1 
4  module grids_io
5 
6  use esmf
7  use netcdf
8  use esmf_logpublicmod
9  use utilities, only : error_handler, netcdf_err
10 
11  implicit none
12 
13  private
14 
15  integer, public, parameter :: n_tiles=6
16  ! mask values for land / ocean mask built from veg type
17  integer, public, parameter :: vtype_water=0, & !< non-land
18  vtype_landice=15
19  ! mask values for soilsnow_mask calculated in the GSI EnKF
20  integer, public, parameter :: mtype_water=0, & !< water
21  mtype_snow=2
22  type, public :: grid_setup_type
23  character(7) :: descriptor
24  character(100) :: fname
25  character(100) :: dir
26  character(15) :: mask_variable(1)
27  character(100) :: fname_mask
28  character(100) :: dir_mask
29  logical :: mask_from_input
30  character(100) :: fname_coord
31  character(100) :: dir_coord
32  integer :: ires
33  integer :: jres
34  end type
35 
36  public :: setup_grid, &
37  read_into_fields, &
38  write_from_fields
39 
40  contains
41 
48 
49  subroutine setup_grid(localpet, npets, grid_setup, mod_grid, timestamp )
50 
51  implicit none
52 
53  ! INTENT IN
54  type(grid_setup_type), intent(in) :: grid_setup
55  integer, intent(in) :: localpet, npets
56  integer, intent(in), optional :: timestamp
57 
58  ! INTENT OUT
59  type(esmf_grid), intent(out) :: mod_grid
60 
61 
62  ! LOCAL
63  type(esmf_field) :: mask_field(1,1)
64  real(esmf_kind_r8), pointer :: ptr_maskvar(:,:)
65  integer(esmf_kind_i4), pointer :: ptr_mask(:,:)
66 
67  integer :: ierr, ncid, tile
68  character(len=128) :: fname_mask
69  character(len=3) :: tstr
70 
71 !--------------------------
72 ! Create grid object, and set up pet distribution
73 
74  select case (grid_setup%descriptor)
75  case ('fv3_rst')
76  call create_grid_fv3(grid_setup%ires, trim(grid_setup%dir_coord), npets, localpet ,mod_grid)
77  case ('gau_inc')
78  call create_grid_gauss(grid_setup, npets, localpet, mod_grid)
79  case default
80  call error_handler("unknown grid_setup%descriptor in setup_grid", 1)
81  end select
82 
83 !--------------------------
84 ! Calculate and add the mask
85 
86  mask_field(1,1) = esmf_fieldcreate(mod_grid, &
87  typekind=esmf_typekind_r8, &
88  staggerloc=esmf_staggerloc_center, &
89  name="input variable for mask", &
90  rc=ierr)
91  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
92  call error_handler("IN FieldCreate, mask_variable", ierr)
93 
94  if (present(timestamp)) then
95  write(tstr,"(I3.3)") timestamp
96  fname_mask = trim(grid_setup%fname_mask)//tstr//".nc"
97  else
98  fname_mask = trim(grid_setup%fname_mask)
99  endif
100 
101  call read_into_fields(localpet, grid_setup%ires, grid_setup%jres, trim(fname_mask), &
102  trim(grid_setup%dir_mask), grid_setup, 1, &
103  grid_setup%mask_variable(1), mask_field(1,1))
104 
105 ! get pointer to mask
106  call esmf_fieldget(mask_field(1,1), &
107  farrayptr=ptr_maskvar, &
108  rc=ierr)
109  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
110  call error_handler("IN FieldGet", ierr)
111 
112 ! create and get pointer to the mask
113  call esmf_gridadditem(mod_grid, &
114  itemflag=esmf_griditem_mask, &
115  staggerloc=esmf_staggerloc_center, &
116  rc=ierr)
117  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
118  call error_handler("in GridAddItem mask", ierr)
119 
120  call esmf_gridgetitem(mod_grid, &
121  itemflag=esmf_griditem_mask, &
122  farrayptr=ptr_mask, &
123  rc=ierr)
124  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
125  call error_handler("in GridGetItem mask", ierr)
126 
127 ! calculate the mask
128  ptr_mask = 1 ! initialize land everywhere
129  select case (trim(grid_setup%mask_variable(1)))
130  case("vegetation_type") ! removing non-land and glaciers using veg class
131  where (nint(ptr_maskvar) == vtype_water ) ptr_mask = 0 ! exclude water
132  where (nint(ptr_maskvar) == vtype_landice ) ptr_mask = 0 ! exclude glaciers
133  case("soilsnow_mask") ! removing snow and non-land using pre-computed mask
134  where (nint(ptr_maskvar) == mtype_water ) ptr_mask = 0 ! exclude non-soil
135  where (nint(ptr_maskvar) == mtype_snow ) ptr_mask = 0 ! exclude snow
136  case default
137  call error_handler("unknown mask_variable", 1)
138  end select
139 
140 ! destroy mask field
141  call esmf_fielddestroy(mask_field(1,1),rc=ierr)
142  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
143  call error_handler("DESTROYING FIELD", ierr)
144 
145  end subroutine setup_grid
146 
157 
158  subroutine read_into_fields(localpet, i_dim, j_dim , fname_read, dir_read, &
159  grid_setup, n_vars, variable_list, fields)
161  implicit none
162 
163  ! INTENT IN
164  integer, intent(in) :: localpet, i_dim, j_dim, n_vars
165  character(*), intent(in) :: fname_read
166  character(*), intent(in) :: dir_read
167  type(grid_setup_type), intent(in) :: grid_setup
168 
169  character(len=15), dimension(n_vars), intent(in) :: variable_list
170 
171  ! INTENT OUT
172  type(esmf_field), dimension(1,n_vars), intent(inout) :: fields
173 
174  ! LOCAL
175  integer :: tt, id_var, ncid, ierr, v, j
176  integer :: n_files
177  character(len=1) :: tchar
178  character(len=500) :: fname
179  real(esmf_kind_r8), allocatable :: array2d(:,:)
180  real(esmf_kind_r8), allocatable :: array_in(:,:,:)
181  real(esmf_kind_r8), allocatable :: temp_array(:,:,:)
182 
183  allocate(array_in(n_vars,i_dim, j_dim))
184  allocate(array2d(i_dim, j_dim))
185 
186  select case (grid_setup%descriptor)
187  case ('fv3_rst')
188  n_files=n_tiles
189  case ('gau_inc')
190  n_files=1
191  case default
192  call error_handler("unknown grid_setup%descriptor in read into fields", 1)
193  end select
194 
195  do tt = 1, n_files
196 
197  ! read from restart
198  if (localpet == 0) then
199 
200  if ( n_files > 1) then
201  write(tchar,'(i1)') tt
202  fname = dir_read//"/"//fname_read//".tile"//tchar//".nc"
203  else
204  fname = dir_read//"/"//fname_read
205  endif
206 
207  print *, 'Reading ', trim(fname)
208 
209  ierr=nf90_open(trim(fname),nf90_nowrite,ncid)
210  call netcdf_err(ierr, 'opening: '//trim(fname) )
211 
212  do v =1, n_vars
213  print *, 'Reading ', trim(variable_list(v))
214  ierr=nf90_inq_varid(ncid, trim(variable_list(v)), id_var)
215  call netcdf_err(ierr, 'reading variable id' )
216 
217  ierr=nf90_get_var(ncid, id_var, array_in(v,:,:))
218  call netcdf_err(ierr, 'reading variable' )
219  enddo
220  ierr = nf90_close(ncid)
221 
222  ! increment files are S->N, ESMF expects N->S
223  if ( grid_setup%descriptor == 'gau_inc') then
224  allocate(temp_array(n_vars,i_dim, j_dim))
225  temp_array = array_in
226  do j=1,j_dim
227  array_in(:,:,j) = temp_array(:,:,j_dim-j+1)
228  enddo
229  deallocate(temp_array)
230  endif
231 
232  endif
233  ! scatter
234  do v =1, n_vars
235  array2d=array_in(v,:,:) ! scatter misbehaves if given indexed 3D array.
236  call esmf_fieldscatter(fields(1,v), array2d, rootpet=0, tile=tt, rc=ierr)
237  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
238  call error_handler("IN FieldScatter", ierr)
239 
240  enddo
241 
242  enddo
243 
244  ! clean up
245  deallocate(array_in)
246  deallocate(array2d)
247 
248  end subroutine read_into_fields
249 
261 
262  subroutine write_from_fields(localpet, i_dim, j_dim , fname_out, dir_out, &
263  n_vars, n_tims, variable_list, fields, add_time_dim)
265  implicit none
266 
267  ! INTENT IN
268  integer, intent(in) :: localpet, i_dim, j_dim, n_vars, n_tims
269  character(*), intent(in) :: fname_out
270  character(*), intent(in) :: dir_out
271  character(15), dimension(n_vars), intent(in) :: variable_list
272  type(esmf_field), dimension(n_tims,n_vars), intent(in) :: fields
273  logical, intent(in) :: add_time_dim
274 
275  ! LOCAL
276  integer :: tt, id_var, ncid, ierr, &
277  id_x, id_y, id_t, v, t
278  character(len=1) :: tchar
279  character(len=500) :: fname
280  real(esmf_kind_r8), allocatable :: array2d(:,:)
281  real(esmf_kind_r8), allocatable :: array_out(:,:,:,:)
282 
283  do v = 1, n_vars
284  if (localpet == 0) print *, 'Writing ', trim(variable_list(v)), ' into field'
285  enddo
286 
287  if (localpet==0) then
288  allocate(array_out(n_vars, i_dim, j_dim, n_tims))
289  allocate(array2d(i_dim, j_dim))
290  else
291  allocate(array_out(0,0,0,0))
292  allocate(array2d(0,0))
293  end if
294 
295  do tt = 1, n_tiles
296 
297  ! fetch data (all PETs)
298  do t =1 , n_tims
299  do v = 1, n_vars
300  call esmf_fieldgather(fields(t,v), array2d, rootpet=0, tile=tt, rc=ierr)
301  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
302  call error_handler("IN FieldGather", ierr)
303  array_out(v,:,:,t) = array2d
304  enddo
305  enddo
306 
307  ! write to netcdf
308  if (localpet == 0) then
309 
310  ! open file, set dimensions
311  write(tchar,'(i1)') tt
312  fname = dir_out//"/"//fname_out//".tile"//tchar//".nc"
313 
314  ierr = nf90_create(trim(fname), nf90_netcdf4, ncid)
315  call netcdf_err(ierr, 'creating file='//trim(fname) )
316 
317  if (add_time_dim) then ! UFS_UTILS expects input with no time dim
318  ! GFS (for IAU) expects a time dimension
319  ! later: update GFS to not expect a time dimension
320  ierr = nf90_def_dim(ncid, 'Time', n_tims, id_t)
321  call netcdf_err(ierr, 'defining taxis dimension' )
322  endif
323 
324  ierr = nf90_def_dim(ncid, 'xaxis_1', i_dim, id_x)
325  call netcdf_err(ierr, 'defining xaxis dimension' )
326 
327  ierr = nf90_def_dim(ncid, 'yaxis_1', j_dim, id_y)
328  call netcdf_err(ierr, 'defining yaxis dimension' )
329 
330 
331  do v=1, n_vars
332 
333  if (add_time_dim) then
334  ! UFS model code to read in the increments is expecting
335  ! dimensions: time, y, x (in the ncdump read out - which reverses fortran indexes)
336  ! need dimensions to be x,y,t below.
337  ierr = nf90_def_var(ncid, trim(variable_list(v)), nf90_double, &
338  (/id_x, id_y, id_t/) , id_var)
339 
340  call netcdf_err(ierr, 'defining '//variable_list(v) )
341  else
342  ierr = nf90_def_var(ncid, trim(variable_list(v)), nf90_double, &
343  (/id_x, id_y/) , id_var)
344  endif
345 
346  call netcdf_err(ierr, 'defining '//variable_list(v) )
347 
348  ierr = nf90_put_var( ncid, id_var, array_out(v,:,:,:) )
349  call netcdf_err(ierr, 'writing '//variable_list(v) )
350 
351  enddo
352 
353  ierr = nf90_close(ncid)
354 
355  endif
356 
357  enddo
358 
359  ! clean up
360  deallocate(array_out)
361 
362  end subroutine write_from_fields
363 
371 
372 
373  subroutine create_grid_fv3(res_atm, dir_fix, npets, localpet, fv3_grid)
375 ! INTENT IN
376  integer, intent(in) :: npets, localpet
377  integer, intent(in) :: res_atm
378  character(*), intent(in) :: dir_fix
379 
380  ! INTENT OUT
381  type(esmf_grid), intent(out) :: fv3_grid
382 
383  integer :: ierr, extra, tile
384  integer :: decomptile(2,n_tiles)
385 
386  character(len=5) :: rchar
387  character(len=200) :: fname
388 
389  if (localpet == 0) print*," creating fv3 grid for ", res_atm
390 
391 ! pet distribution
392  extra = npets / n_tiles
393  do tile = 1, n_tiles
394  decomptile(:,tile)=(/1,extra/)
395  enddo
396 
397  ! mosaic file
398  write(rchar,'(i5)') res_atm
399  fname = trim(dir_fix)//"/C"//trim(adjustl(rchar))// "_mosaic.nc"
400 
401 ! create the grid
402  fv3_grid = esmf_gridcreatemosaic(filename=trim(fname), &
403  regdecompptile=decomptile, &
404  staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
405  esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
406  indexflag=esmf_index_global, &
407  tilefilepath=trim(dir_fix), &
408  rc=ierr)
409  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
410  call error_handler("IN GridCreateMosaic", ierr)
411 
412  end subroutine create_grid_fv3
413 
420 
421  subroutine create_grid_gauss(grid_setup, npets, localpet, gauss_grid)
423  ! INTENT IN
424  type(grid_setup_type), intent(in) :: grid_setup
425  integer, intent(in) :: npets, localpet
426 
427  ! INTENT OUT
428  type(esmf_grid) :: gauss_grid
429 
430  integer :: ierr, fac
431  character(len=200) :: fname
432 
433  fname = trim(grid_setup%dir_coord)//trim(grid_setup%fname_coord)
434 
435  if (localpet == 0) print*," creating gauss grid for ", trim(fname)
436 
437  fac = npets / n_tiles
438  gauss_grid = esmf_gridcreate(filename=trim(fname), &
439  fileformat=esmf_fileformat_scrip, &
440  regdecomp=(/n_tiles,fac/), addcornerstagger=.true., rc=ierr)
441  if(esmf_logfounderror(rctocheck=ierr,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
442  call error_handler("IN Gauss GridCreate", ierr)
443 
444  end subroutine create_grid_gauss
445 
446 end module grids_io