FV3-JEDI
fv3jedi_geom_mod.f90
Go to the documentation of this file.
1 ! (C) Copyright 2017-2020 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 
6 !> Fortran module handling geometry for the FV3 model
7 
8 module fv3jedi_geom_mod
9 
10 use netcdf
11 use mpi
12 use string_f_c_mod
13 
14 ! atlas uses
15 use atlas_module, only: atlas_field, atlas_fieldset, atlas_real, atlas_functionspace_pointcloud
16 
17 ! fckit uses
18 use fckit_mpi_module, only: fckit_mpi_comm
19 use fckit_configuration_module, only: fckit_configuration
20 
21 ! fms uses
22 use fms_io_mod, only: fms_io_init, fms_io_exit
23 use fms_mod, only: fms_init
24 use mpp_mod, only: mpp_init, mpp_exit, mpp_pe, mpp_npes, mpp_error, fatal, note
25 use mpp_domains_mod, only: domain2d, mpp_deallocate_domain, mpp_domains_init, &
26  mpp_define_layout, mpp_define_mosaic, mpp_define_io_domain, &
27  mpp_domains_exit, mpp_domains_set_stack_size
28 
29 ! fv3 uses
30 use fv_arrays_mod, only: fv_atmos_type, deallocate_fv_atmos_type
31 
32 ! fv3jedi uses
33 use fields_metadata_mod, only: fields_metadata, field_metadata
34 use fv3jedi_constants_mod, only: ps, rad2deg
35 use fv3jedi_kinds_mod, only: kind_real
36 use fv3jedi_netcdf_utils_mod, only: nccheck
37 use fv_init_mod, only: fv_init
38 
39 implicit none
40 private
41 public :: fv3jedi_geom, getverticalcoordlogp, initialize
42 
43 ! --------------------------------------------------------------------------------------------------
44 
45 !> Fortran derived type to hold geometry data for the FV3JEDI model
46 type :: fv3jedi_geom
47  integer :: isd, ied, jsd, jed !data domain
48  integer :: isc, iec, jsc, jec !compute domain
49  integer :: npx,npy,npz,ngrid !x/y/z-dir grid edge points per tile
50  integer :: layout(2), io_layout(2) !Processor layouts
51  integer :: ntile, ntiles !Tile number and total
52  real(kind=kind_real) :: ptop !Pressure at top of domain
53  type(domain2d) :: domain_fix !MPP domain
54  type(domain2d), pointer :: domain !MPP domain
55  character(len=10) :: interp_method !Interpolation type
56  real(kind=kind_real), allocatable, dimension(:) :: ak, bk !Model level coefficients
57  real(kind=kind_real), allocatable, dimension(:,:) :: grid_lon, grid_lat !Lat/lon centers
58  real(kind=kind_real), allocatable, dimension(:,:) :: egrid_lon, egrid_lat !Lat/lon edges
59  real(kind=kind_real), allocatable, dimension(:) :: lon_us, lat_us !Lat/lon centers unstructured
60  real(kind=kind_real), allocatable, dimension(:,:) :: area !Grid area
61  real(kind=kind_real), allocatable, dimension(:,:) :: dx, dy !dx/dy at edges
62  real(kind=kind_real), allocatable, dimension(:,:) :: dxc, dyc !dx/dy c grid
63  real(kind=kind_real), allocatable, dimension(:,:,:) :: grid, vlon, vlat
64  real(kind=kind_real), allocatable, dimension(:) :: edge_vect_n, edge_vect_e
65  real(kind=kind_real), allocatable, dimension(:) :: edge_vect_s, edge_vect_w
66  real(kind=kind_real), allocatable, dimension(:,:,:,:) :: es, ew
67  real(kind=kind_real), allocatable, dimension(:,:) :: a11, a12, a21, a22
68  type(fckit_mpi_comm) :: f_comm
69  type(fields_metadata) :: fields
70  ! For D to (A to) C grid
71  real(kind=kind_real), allocatable, dimension(:,:) :: rarea
72  real(kind=kind_real), allocatable, dimension(:,:,:) :: sin_sg
73  real(kind=kind_real), allocatable, dimension(:,:) :: cosa_u
74  real(kind=kind_real), allocatable, dimension(:,:) :: cosa_v
75  real(kind=kind_real), allocatable, dimension(:,:) :: cosa_s
76  real(kind=kind_real), allocatable, dimension(:,:) :: rsin_u
77  real(kind=kind_real), allocatable, dimension(:,:) :: rsin_v
78  real(kind=kind_real), allocatable, dimension(:,:) :: rsin2
79  real(kind=kind_real), allocatable, dimension(:,:) :: dxa, dya
80  logical :: ne_corner, se_corner, sw_corner, nw_corner
81  logical :: nested = .false.
82  integer :: grid_type = 0
83  logical :: dord4 = .true.
84  type(atlas_functionspace_pointcloud) :: afunctionspace
85  contains
86  procedure, public :: create
87  procedure, public :: clone
88  procedure, public :: delete
89  procedure, public :: set_atlas_lonlat
90  procedure, public :: fill_atlas_fieldset
91 end type fv3jedi_geom
92 
93 ! --------------------------------------------------------------------------------------------------
94 
95 contains
96 
97 ! --------------------------------------------------------------------------------------------------
98 
99 subroutine initialize(conf, comm)
100 
101 type(fckit_configuration), intent(in) :: conf
102 type(fckit_mpi_comm), intent(in) :: comm
103 
104 
105 integer :: stackmax = 4000000
106 
107 call mpp_init(localcomm=comm%communicator())
108 call mpp_domains_init
109 call fms_io_init
110 call fms_init()
111 
112 if (conf%has("stackmax")) call conf%get_or_die("stackmax",stackmax)
113 
114 call mpp_domains_set_stack_size(stackmax)
115 
116 end subroutine initialize
117 
118 ! --------------------------------------------------------------------------------------------------
119 
120 subroutine create(self, conf, comm, fields)
121 
122 !Arguments
123 class(fv3jedi_geom), target, intent(inout) :: self
124 type(fckit_configuration), intent(in) :: conf
125 type(fckit_mpi_comm), intent(in) :: comm
126 type(fields_metadata), intent(in) :: fields
127 
128 !Locals
129 character(len=256) :: file_akbk
130 type(fv_atmos_type), allocatable :: Atm(:)
131 logical, allocatable :: grids_on_this_pe(:)
132 integer :: i, j, jj, gtile
133 integer :: p_split = 1
134 integer :: ncstat, ncid, akvarid, bkvarid, readdim, dcount
135 integer, dimension(nf90_max_var_dims) :: dimids, dimlens
136 
137 character(len=:), allocatable :: str
138 logical :: do_write_geom = .false.
139 
140 ! Add the communicator to the geometry
141 ! ------------------------------------
142 self%f_comm = comm
143 
144 ! Interpolation type
145 ! ------------------
146 self%interp_method = 'barycent'
147 if (conf%has("interpolation method")) then
148  call conf%get_or_die("interpolation method",str)
149  self%interp_method = str
150  deallocate(str)
151 endif
152 
153 ! Local pointer to field meta data
154 ! --------------------------------
155 self%fields = fields
156 
157 !Intialize using the model setup routine
158 ! --------------------------------------
159 call fv_init(atm, 300.0_kind_real, grids_on_this_pe, p_split, gtile)
160 
161 self%isd = atm(1)%bd%isd
162 self%ied = atm(1)%bd%ied
163 self%jsd = atm(1)%bd%jsd
164 self%jed = atm(1)%bd%jed
165 self%isc = atm(1)%bd%isc
166 self%iec = atm(1)%bd%iec
167 self%jsc = atm(1)%bd%jsc
168 self%jec = atm(1)%bd%jec
169 
170 self%ntile = gtile
171 self%ntiles = atm(1)%flagstruct%ntiles
172 
173 self%npx = atm(1)%npx
174 self%npy = atm(1)%npy
175 self%npz = atm(1)%npz
176 
177 self%layout(1) = atm(1)%layout(1)
178 self%layout(2) = atm(1)%layout(2)
179 self%io_layout(1) = atm(1)%io_layout(1)
180 self%io_layout(2) = atm(1)%io_layout(2)
181 
182 !Allocatable arrays
183 allocate(self%ak(self%npz+1) )
184 allocate(self%bk(self%npz+1) )
185 
186 allocate(self%grid_lon (self%isd :self%ied, self%jsd :self%jed ))
187 allocate(self%grid_lat (self%isd :self%ied, self%jsd :self%jed ))
188 allocate(self%egrid_lon (self%isd :self%ied+1,self%jsd :self%jed+1))
189 allocate(self%egrid_lat (self%isd :self%ied+1,self%jsd :self%jed+1))
190 allocate(self%area (self%isd :self%ied, self%jsd :self%jed ))
191 allocate(self%dx (self%isd :self%ied ,self%jsd :self%jed+1))
192 allocate(self%dy (self%isd :self%ied+1,self%jsd :self%jed ))
193 allocate(self%dxc (self%isd :self%ied+1,self%jsd :self%jed ))
194 allocate(self%dyc (self%isd :self%ied ,self%jsd :self%jed+1))
195 
196 allocate(self%grid (self%isd :self%ied+1,self%jsd :self%jed+1,2))
197 allocate(self%vlon (self%isc-2:self%iec+2,self%jsc-2:self%jec+2,3))
198 allocate(self%vlat (self%isc-2:self%iec+2,self%jsc-2:self%jec+2,3))
199 
200 allocate(self%edge_vect_n(self%isd:self%ied))
201 allocate(self%edge_vect_e(self%jsd:self%jed))
202 allocate(self%edge_vect_s(self%isd:self%ied))
203 allocate(self%edge_vect_w(self%jsd:self%jed))
204 
205 allocate(self%es(3,self%isd:self%ied ,self%jsd:self%jed+1,2))
206 allocate(self%ew(3,self%isd:self%ied+1,self%jsd:self%jed, 2))
207 
208 allocate(self%a11(self%isc-1:self%iec+1,self%jsc-1:self%jec+1) )
209 allocate(self%a12(self%isc-1:self%iec+1,self%jsc-1:self%jec+1) )
210 allocate(self%a21(self%isc-1:self%iec+1,self%jsc-1:self%jec+1) )
211 allocate(self%a22(self%isc-1:self%iec+1,self%jsc-1:self%jec+1) )
212 
213 allocate(self%rarea (self%isd:self%ied ,self%jsd:self%jed ))
214 allocate(self%sin_sg(self%isd:self%ied ,self%jsd:self%jed ,9))
215 allocate(self%cosa_u(self%isd:self%ied+1,self%jsd:self%jed ))
216 allocate(self%cosa_v(self%isd:self%ied ,self%jsd:self%jed+1))
217 allocate(self%cosa_s(self%isd:self%ied ,self%jsd:self%jed ))
218 allocate(self%rsin_u(self%isd:self%ied+1,self%jsd:self%jed ))
219 allocate(self%rsin_v(self%isd:self%ied ,self%jsd:self%jed+1))
220 allocate(self%rsin2 (self%isd:self%ied ,self%jsd:self%jed ))
221 allocate(self%dxa (self%isd:self%ied ,self%jsd:self%jed ))
222 allocate(self%dya (self%isd:self%ied ,self%jsd:self%jed ))
223 
224 ! ak and bk hybrid coordinate coefficients
225 ! ----------------------------------------
226 if (self%npz > 1) then
227 
228  ! Set path/filename for ak and bk file
229  call conf%get_or_die("akbk",str)
230  file_akbk = str
231 
232  !Open file
233  call nccheck ( nf90_open(file_akbk, nf90_nowrite, ncid), "nf90_open "//file_akbk )
234 
235  !Search for ak in the file
236  ncstat = nf90_inq_varid(ncid, "ak", akvarid)
237  if(ncstat /= nf90_noerr) call abor1_ftn("Failed to find ak in file "//file_akbk)
238 
239  !Search for bk in the file
240  ncstat = nf90_inq_varid(ncid, "bk", bkvarid)
241  if(ncstat /= nf90_noerr) call abor1_ftn("Failed to find bk in file "//file_akbk)
242 
243  ! Check that dimension of ak/bk in the file match vertical levels of model
244  dimids = 0
245  call nccheck ( nf90_inquire_variable(ncid, akvarid, dimids = dimids), "nf90_inq_var ak" )
246  readdim = -1
247  dcount = 0
248  do i = 1,nf90_max_var_dims
249  if (dimids(i) > 0) then
250  call nccheck( nf90_inquire_dimension(ncid, dimids(i), len = dimlens(i)), &
251  "nf90_inquire_dimension" )
252  if (dimlens(i) == self%npz+1) then
253  readdim = i
254  endif
255  dcount = dcount + 1
256  endif
257  enddo
258  if (readdim == -1) call abor1_ftn("ak/bk in file does not match dimension of npz from input.nml")
259 
260  !Read ak and bk from the file
261  call nccheck( nf90_get_var(ncid, akvarid, self%ak), "fv3jedi_geom, nf90_get_var ak" )
262  call nccheck( nf90_get_var(ncid, bkvarid, self%bk), "fv3jedi_geom, nf90_get_var bk" )
263 else
264  self%ak = 0.0_kind_real
265  self%bk = 0.0_kind_real
266 endif
267 
268 ! Arrays from the Atm Structure
269 ! --------------------------------
270 
271 self%grid_lon = real(atm(1)%gridstruct%agrid_64(:,:,1),kind_real)
272 self%grid_lat = real(atm(1)%gridstruct%agrid_64(:,:,2),kind_real)
273 self%egrid_lon = real(atm(1)%gridstruct%grid_64(:,:,1),kind_real)
274 self%egrid_lat = real(atm(1)%gridstruct%grid_64(:,:,2),kind_real)
275 self%area = real(atm(1)%gridstruct%area_64,kind_real)
276 self%dx = real(atm(1)%gridstruct%dx ,kind_real)
277 self%dy = real(atm(1)%gridstruct%dy ,kind_real)
278 self%dxc = real(atm(1)%gridstruct%dxc,kind_real)
279 self%dyc = real(atm(1)%gridstruct%dyc,kind_real)
280 
281 self%grid = real(atm(1)%gridstruct%grid,kind_real)
282 self%vlon = real(atm(1)%gridstruct%vlon,kind_real)
283 self%vlat = real(atm(1)%gridstruct%vlat,kind_real)
284 
285 self%edge_vect_n = real(atm(1)%gridstruct%edge_vect_n,kind_real)
286 self%edge_vect_e = real(atm(1)%gridstruct%edge_vect_e,kind_real)
287 self%edge_vect_s = real(atm(1)%gridstruct%edge_vect_s,kind_real)
288 self%edge_vect_w = real(atm(1)%gridstruct%edge_vect_w,kind_real)
289 
290 self%es = real(atm(1)%gridstruct%es,kind_real)
291 self%ew = real(atm(1)%gridstruct%ew,kind_real)
292 
293 self%a11 = real(atm(1)%gridstruct%a11,kind_real)
294 self%a12 = real(atm(1)%gridstruct%a12,kind_real)
295 self%a21 = real(atm(1)%gridstruct%a21,kind_real)
296 self%a22 = real(atm(1)%gridstruct%a22,kind_real)
297 
298 self%rarea = real(atm(1)%gridstruct%rarea ,kind_real)
299 self%sin_sg = real(atm(1)%gridstruct%sin_sg,kind_real)
300 self%cosa_u = real(atm(1)%gridstruct%cosa_u,kind_real)
301 self%cosa_v = real(atm(1)%gridstruct%cosa_v,kind_real)
302 self%cosa_s = real(atm(1)%gridstruct%cosa_s,kind_real)
303 self%rsin_u = real(atm(1)%gridstruct%rsin_u,kind_real)
304 self%rsin_v = real(atm(1)%gridstruct%rsin_v,kind_real)
305 self%rsin2 = real(atm(1)%gridstruct%rsin2 ,kind_real)
306 self%dxa = real(atm(1)%gridstruct%dxa ,kind_real)
307 self%dya = real(atm(1)%gridstruct%dya ,kind_real)
308 self%ne_corner = atm(1)%gridstruct%ne_corner
309 self%se_corner = atm(1)%gridstruct%se_corner
310 self%sw_corner = atm(1)%gridstruct%sw_corner
311 self%nw_corner = atm(1)%gridstruct%nw_corner
312 
313 !Unstructured lat/lon
314 self%ngrid = (self%iec-self%isc+1)*(self%jec-self%jsc+1)
315 allocate(self%lat_us(self%ngrid))
316 allocate(self%lon_us(self%ngrid))
317 
318 jj = 0
319 do j = self%jsc,self%jec
320  do i = self%isc,self%iec
321  jj = jj + 1
322  self%lat_us(jj) = self%grid_lat(i,j)
323  self%lon_us(jj) = self%grid_lon(i,j)
324  enddo
325 enddo
326 
327 !Set Ptop
328 self%ptop = self%ak(1)
329 
330 !Done with the Atm stucture here
331 call deallocate_fv_atmos_type(atm(1))
332 deallocate(atm)
333 deallocate(grids_on_this_pe)
334 
335 !Resetup domain to avoid risk of copied pointers
336 call setup_domain( self%domain_fix, self%npx-1, self%npy-1, &
337  self%ntiles, self%layout, self%io_layout, 3)
338 
339 self%domain => self%domain_fix
340 
341 ! Optionally write the geometry to file
342 ! -------------------------------------
343 if (conf%has("do_write_geom")) then
344  call conf%get_or_die("do_write_geom",do_write_geom)
345 endif
346 
347 if (do_write_geom) then
348  call write_geom(self)
349 endif
350 
351 end subroutine create
352 
353 ! --------------------------------------------------------------------------------------------------
354 
355 subroutine clone(self, other, fields)
356 
357 class(fv3jedi_geom), intent(inout) :: self
358 type(fv3jedi_geom), target, intent(in) :: other
359 type(fields_metadata), intent(in) :: fields
360 
361 allocate(self%ak(other%npz+1) )
362 allocate(self%bk(other%npz+1) )
363 
364 allocate(self%grid_lon (other%isd :other%ied, other%jsd :other%jed ))
365 allocate(self%grid_lat (other%isd :other%ied, other%jsd :other%jed ))
366 allocate(self%egrid_lon (other%isd :other%ied+1,other%jsd :other%jed+1))
367 allocate(self%egrid_lat (other%isd :other%ied+1,other%jsd :other%jed+1))
368 allocate(self%area (other%isd :other%ied, other%jsd :other%jed ))
369 allocate(self%dx (other%isd :other%ied ,other%jsd :other%jed+1))
370 allocate(self%dy (other%isd :other%ied+1,other%jsd :other%jed ))
371 allocate(self%dxc (other%isd :other%ied+1,other%jsd :other%jed ))
372 allocate(self%dyc (other%isd :other%ied ,other%jsd :other%jed+1))
373 
374 allocate(self%grid (other%isd :other%ied+1,other%jsd :other%jed+1,2))
375 allocate(self%vlon (other%isc-2:other%iec+2,other%jsc-2:other%jec+2,3))
376 allocate(self%vlat (other%isc-2:other%iec+2,other%jsc-2:other%jec+2,3))
377 
378 allocate(self%edge_vect_n(other%isd:other%ied))
379 allocate(self%edge_vect_e(other%jsd:other%jed))
380 allocate(self%edge_vect_s(other%isd:other%ied))
381 allocate(self%edge_vect_w(other%jsd:other%jed))
382 
383 allocate(self%es(3,other%isd:other%ied ,other%jsd:other%jed+1,2))
384 allocate(self%ew(3,other%isd:other%ied+1,other%jsd:other%jed, 2))
385 
386 allocate(self%a11(other%isc-1:other%iec+1,other%jsc-1:other%jec+1) )
387 allocate(self%a12(other%isc-1:other%iec+1,other%jsc-1:other%jec+1) )
388 allocate(self%a21(other%isc-1:other%iec+1,other%jsc-1:other%jec+1) )
389 allocate(self%a22(other%isc-1:other%iec+1,other%jsc-1:other%jec+1) )
390 
391 allocate(self%rarea (other%isd:other%ied ,other%jsd:other%jed ))
392 allocate(self%sin_sg(other%isd:other%ied ,other%jsd:other%jed ,9))
393 allocate(self%cosa_u(other%isd:other%ied+1,other%jsd:other%jed ))
394 allocate(self%cosa_v(other%isd:other%ied ,other%jsd:other%jed+1))
395 allocate(self%cosa_s(other%isd:other%ied ,other%jsd:other%jed ))
396 allocate(self%rsin_u(other%isd:other%ied+1,other%jsd:other%jed ))
397 allocate(self%rsin_v(other%isd:other%ied ,other%jsd:other%jed+1))
398 allocate(self%rsin2 (other%isd:other%ied ,other%jsd:other%jed ))
399 allocate(self%dxa (other%isd:other%ied ,other%jsd:other%jed ))
400 allocate(self%dya (other%isd:other%ied ,other%jsd:other%jed ))
401 
402 allocate(self%lat_us(other%ngrid))
403 allocate(self%lon_us(other%ngrid))
404 
405 self%npx = other%npx
406 self%npy = other%npy
407 self%npz = other%npz
408 self%ngrid = other%ngrid
409 self%layout = other%layout
410 self%io_layout = other%io_layout
411 self%isc = other%isc
412 self%isd = other%isd
413 self%iec = other%iec
414 self%ied = other%ied
415 self%jsc = other%jsc
416 self%jsd = other%jsd
417 self%jec = other%jec
418 self%jed = other%jed
419 self%ntile = other%ntile
420 self%ntiles = other%ntiles
421 self%ptop = other%ptop
422 self%ak = other%ak
423 self%bk = other%bk
424 self%grid_lon = other%grid_lon
425 self%grid_lat = other%grid_lat
426 self%egrid_lon = other%egrid_lon
427 self%egrid_lat = other%egrid_lat
428 self%area = other%area
429 self%dx = other%dx
430 self%dy = other%dy
431 self%dxc = other%dxc
432 self%dyc = other%dyc
433 self%grid = other%grid
434 self%vlon = other%vlon
435 self%vlat = other%vlat
436 self%edge_vect_n = other%edge_vect_n
437 self%edge_vect_e = other%edge_vect_e
438 self%edge_vect_s = other%edge_vect_s
439 self%edge_vect_w = other%edge_vect_w
440 self%es = other%es
441 self%ew = other%ew
442 self%a11 = other%a11
443 self%a12 = other%a12
444 self%a21 = other%a21
445 self%a22 = other%a22
446 self%f_comm = other%f_comm
447 self%interp_method = other%interp_method
448 
449 self%rarea = other%rarea
450 self%sin_sg = other%sin_sg
451 self%cosa_u = other%cosa_u
452 self%cosa_v = other%cosa_v
453 self%cosa_s = other%cosa_s
454 self%rsin_u = other%rsin_u
455 self%rsin_v = other%rsin_v
456 self%rsin2 = other%rsin2
457 self%dxa = other%dxa
458 self%dya = other%dya
459 self%ne_corner = other%ne_corner
460 self%se_corner = other%se_corner
461 self%sw_corner = other%sw_corner
462 self%nw_corner = other%nw_corner
463 
464 self%domain => other%domain
465 
466 self%afunctionspace = atlas_functionspace_pointcloud(other%afunctionspace%c_ptr())
467 
468 self%fields = fields
469 
470 self%lat_us = other%lat_us
471 self%lon_us = other%lon_us
472 
473 end subroutine clone
474 
475 ! --------------------------------------------------------------------------------------------------
476 
477 subroutine delete(self)
478 
479 class(fv3jedi_geom), intent(inout) :: self
480 
481 ! Deallocate
482 deallocate(self%ak)
483 deallocate(self%bk)
484 deallocate(self%grid_lon)
485 deallocate(self%grid_lat)
486 deallocate(self%egrid_lon)
487 deallocate(self%egrid_lat)
488 deallocate(self%area)
489 deallocate(self%dx)
490 deallocate(self%dy)
491 deallocate(self%dxc)
492 deallocate(self%dyc)
493 deallocate(self%grid)
494 deallocate(self%vlon)
495 deallocate(self%vlat)
496 deallocate(self%edge_vect_n)
497 deallocate(self%edge_vect_e)
498 deallocate(self%edge_vect_s)
499 deallocate(self%edge_vect_w)
500 deallocate(self%es)
501 deallocate(self%ew)
502 deallocate(self%a11)
503 deallocate(self%a12)
504 deallocate(self%a21)
505 deallocate(self%a22)
506 
507 deallocate(self%rarea)
508 deallocate(self%sin_sg)
509 deallocate(self%cosa_u)
510 deallocate(self%cosa_v)
511 deallocate(self%cosa_s)
512 deallocate(self%rsin_u)
513 deallocate(self%rsin_v)
514 deallocate(self%rsin2 )
515 deallocate(self%dxa )
516 deallocate(self%dya )
517 
518 deallocate(self%lat_us)
519 deallocate(self%lon_us)
520 
521 ! Required memory leak, since copying this causes problems
522 !call mpp_deallocate_domain(self%domain_fix)
523 
524 call self%afunctionspace%final()
525 
526 ! Could finalize the fms routines. Possibly needs to be done only when key = 0
527 !call fms_io_exit
528 !call mpp_domains_exit
529 !call mpp_exit
530 
531 end subroutine delete
532 
533 ! --------------------------------------------------------------------------------------------------
534 
535 subroutine set_atlas_lonlat(self, afieldset)
536 
537 !Arguments
538 class(fv3jedi_geom), intent(inout) :: self
539 type(atlas_fieldset), intent(inout) :: afieldset
540 
541 !Locals
542 real(kind_real), pointer :: real_ptr(:,:)
543 type(atlas_field) :: afield
544 
545 ! Create lon/lat field
546 afield = atlas_field(name="lonlat", kind=atlas_real(kind_real), shape=(/2,(self%iec-self%isc+1)*(self%jec-self%jsc+1)/))
547 call afield%data(real_ptr)
548 real_ptr(1,:) = rad2deg*pack(self%grid_lon(self%isc:self%iec,self%jsc:self%jec),.true.)
549 real_ptr(2,:) = rad2deg*pack(self%grid_lat(self%isc:self%iec,self%jsc:self%jec),.true.)
550 call afieldset%add(afield)
551 
552 end subroutine set_atlas_lonlat
553 
554 ! --------------------------------------------------------------------------------------------------
555 
556 subroutine fill_atlas_fieldset(self, afieldset)
557 
558 !Arguments
559 class(fv3jedi_geom), intent(inout) :: self
560 type(atlas_fieldset), intent(inout) :: afieldset
561 
562 !Locals
563 integer :: jl
564 real(kind=kind_real) :: sigmaup, sigmadn
565 real(kind=kind_real), pointer :: real_ptr_1(:), real_ptr_2(:,:)
566 type(atlas_field) :: afield
567 
568 ! Add area
569 afield = self%afunctionspace%create_field(name='area', kind=atlas_real(kind_real), levels=0)
570 call afield%data(real_ptr_1)
571 real_ptr_1 = pack(self%area(self%isc:self%iec,self%jsc:self%jec),.true.)
572 call afieldset%add(afield)
573 call afield%final()
574 
575 ! Add vertical unit
576 afield = self%afunctionspace%create_field(name='vunit', kind=atlas_real(kind_real), levels=self%npz)
577 call afield%data(real_ptr_2)
578 do jl=1,self%npz
579  sigmaup = self%ak(jl+1)/ps+self%bk(jl+1) ! si are now sigmas
580  sigmadn = self%ak(jl )/ps+self%bk(jl )
581  real_ptr_2(jl,:) = 0.5*(sigmaup+sigmadn) ! 'fake' sigma coordinates
582 enddo
583 call afieldset%add(afield)
584 call afield%final()
585 
586 end subroutine fill_atlas_fieldset
587 
588 ! --------------------------------------------------------------------------------------------------
589 
590 subroutine setup_domain(domain, nx, ny, ntiles, layout_in, io_layout, halo)
591 
592  type(domain2d), intent(inout) :: domain
593  integer, intent(in) :: nx, ny, ntiles
594  integer, intent(in) :: layout_in(:), io_layout(:)
595  integer, intent(in) :: halo
596 
597  integer :: pe, npes, npes_per_tile, tile
598  integer :: num_contact, num_alloc
599  integer :: n, layout(2)
600  integer, allocatable, dimension(:,:) :: global_indices, layout2D
601  integer, allocatable, dimension(:) :: pe_start, pe_end
602  integer, allocatable, dimension(:) :: tile1, tile2
603  integer, allocatable, dimension(:) :: istart1, iend1, jstart1, jend1
604  integer, allocatable, dimension(:) :: istart2, iend2, jstart2, jend2
605  integer, allocatable :: tile_id(:)
606  logical :: is_symmetry
607 
608  pe = mpp_pe()
609  npes = mpp_npes()
610 
611  if (mod(npes,ntiles) /= 0) then
612  call mpp_error(note, "setup_domain: npes can not be divided by ntiles")
613  return
614  endif
615  npes_per_tile = npes/ntiles
616  tile = pe/npes_per_tile + 1
617 
618  if (layout_in(1)*layout_in(2) == npes_per_tile) then
619  layout = layout_in
620  else
621  call mpp_define_layout( (/1,nx,1,ny/), npes_per_tile, layout )
622  endif
623 
624  if (io_layout(1) <1 .or. io_layout(2) <1) call mpp_error(fatal, &
625  "setup_domain: both elements of variable io_layout must be positive integer")
626  if (mod(layout(1), io_layout(1)) /= 0 ) call mpp_error(fatal, &
627  "setup_domain: layout(1) must be divided by io_layout(1)")
628  if (mod(layout(2), io_layout(2)) /= 0 ) call mpp_error(fatal, &
629  "setup_domain: layout(2) must be divided by io_layout(2)")
630 
631  allocate(global_indices(4,ntiles), layout2d(2,ntiles), pe_start(ntiles), pe_end(ntiles) )
632 
633  ! select case based off of 1 or 6 tiles
634  select case(ntiles)
635  case ( 1 ) ! FV3-SAR
636  num_contact = 0
637  case ( 6 ) ! FV3 global
638  num_contact = 12
639  case default
640  call mpp_error(fatal, "setup_domain: ntiles != 1 or 6")
641  end select
642 
643  do n = 1, ntiles
644  global_indices(:,n) = (/1,nx,1,ny/)
645  layout2d(:,n) = layout
646  pe_start(n) = (n-1)*npes_per_tile
647  pe_end(n) = n*npes_per_tile-1
648  enddo
649  num_alloc = max(1, num_contact)
650  ! this code copied from domain_decomp in fv_mp_mod.f90
651  allocate(tile1(num_alloc), tile2(num_alloc) )
652  allocate(tile_id(ntiles))
653  allocate(istart1(num_alloc), iend1(num_alloc), jstart1(num_alloc), jend1(num_alloc) )
654  allocate(istart2(num_alloc), iend2(num_alloc), jstart2(num_alloc), jend2(num_alloc) )
655  ! select case based off of 1 or 6 tiles
656  select case(ntiles)
657  case ( 1 ) ! FV3-SAR
658  ! No contacts, do nothing
659  case ( 6 ) ! FV3 global
660  !--- Contact line 1, between tile 1 (EAST) and tile 2 (WEST)
661  tile1(1) = 1; tile2(1) = 2
662  istart1(1) = nx; iend1(1) = nx; jstart1(1) = 1; jend1(1) = ny
663  istart2(1) = 1; iend2(1) = 1; jstart2(1) = 1; jend2(1) = ny
664  !--- Contact line 2, between tile 1 (NORTH) and tile 3 (WEST)
665  tile1(2) = 1; tile2(2) = 3
666  istart1(2) = 1; iend1(2) = nx; jstart1(2) = ny; jend1(2) = ny
667  istart2(2) = 1; iend2(2) = 1; jstart2(2) = ny; jend2(2) = 1
668  !--- Contact line 3, between tile 1 (WEST) and tile 5 (NORTH)
669  tile1(3) = 1; tile2(3) = 5
670  istart1(3) = 1; iend1(3) = 1; jstart1(3) = 1; jend1(3) = ny
671  istart2(3) = nx; iend2(3) = 1; jstart2(3) = ny; jend2(3) = ny
672  !--- Contact line 4, between tile 1 (SOUTH) and tile 6 (NORTH)
673  tile1(4) = 1; tile2(4) = 6
674  istart1(4) = 1; iend1(4) = nx; jstart1(4) = 1; jend1(4) = 1
675  istart2(4) = 1; iend2(4) = nx; jstart2(4) = ny; jend2(4) = ny
676  !--- Contact line 5, between tile 2 (NORTH) and tile 3 (SOUTH)
677  tile1(5) = 2; tile2(5) = 3
678  istart1(5) = 1; iend1(5) = nx; jstart1(5) = ny; jend1(5) = ny
679  istart2(5) = 1; iend2(5) = nx; jstart2(5) = 1; jend2(5) = 1
680  !--- Contact line 6, between tile 2 (EAST) and tile 4 (SOUTH)
681  tile1(6) = 2; tile2(6) = 4
682  istart1(6) = nx; iend1(6) = nx; jstart1(6) = 1; jend1(6) = ny
683  istart2(6) = nx; iend2(6) = 1; jstart2(6) = 1; jend2(6) = 1
684  !--- Contact line 7, between tile 2 (SOUTH) and tile 6 (EAST)
685  tile1(7) = 2; tile2(7) = 6
686  istart1(7) = 1; iend1(7) = nx; jstart1(7) = 1; jend1(7) = 1
687  istart2(7) = nx; iend2(7) = nx; jstart2(7) = ny; jend2(7) = 1
688  !--- Contact line 8, between tile 3 (EAST) and tile 4 (WEST)
689  tile1(8) = 3; tile2(8) = 4
690  istart1(8) = nx; iend1(8) = nx; jstart1(8) = 1; jend1(8) = ny
691  istart2(8) = 1; iend2(8) = 1; jstart2(8) = 1; jend2(8) = ny
692  !--- Contact line 9, between tile 3 (NORTH) and tile 5 (WEST)
693  tile1(9) = 3; tile2(9) = 5
694  istart1(9) = 1; iend1(9) = nx; jstart1(9) = ny; jend1(9) = ny
695  istart2(9) = 1; iend2(9) = 1; jstart2(9) = ny; jend2(9) = 1
696  !--- Contact line 10, between tile 4 (NORTH) and tile 5 (SOUTH)
697  tile1(10) = 4; tile2(10) = 5
698  istart1(10) = 1; iend1(10) = nx; jstart1(10) = ny; jend1(10) = ny
699  istart2(10) = 1; iend2(10) = nx; jstart2(10) = 1; jend2(10) = 1
700  !--- Contact line 11, between tile 4 (EAST) and tile 6 (SOUTH)
701  tile1(11) = 4; tile2(11) = 6
702  istart1(11) = nx; iend1(11) = nx; jstart1(11) = 1; jend1(11) = ny
703  istart2(11) = nx; iend2(11) = 1; jstart2(11) = 1; jend2(11) = 1
704  !--- Contact line 12, between tile 5 (EAST) and tile 6 (WEST)
705  tile1(12) = 5; tile2(12) = 6
706  istart1(12) = nx; iend1(12) = nx; jstart1(12) = 1; jend1(12) = ny
707  istart2(12) = 1; iend2(12) = 1; jstart2(12) = 1; jend2(12) = ny
708  end select
709  is_symmetry = .true.
710  do n = 1, ntiles
711  tile_id(n) = n
712  enddo
713 
714  call mpp_define_mosaic(global_indices, layout2d, domain, ntiles, num_contact, tile1, tile2, &
715  istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2, &
716  pe_start, pe_end, whalo=halo, ehalo=halo, shalo=halo, nhalo=halo, &
717  symmetry=is_symmetry, tile_id=tile_id, &
718  name='cubic_grid')
719 
720  if (io_layout(1) /= 1 .or. io_layout(2) /= 1) call mpp_define_io_domain(domain, io_layout)
721 
722  deallocate(pe_start, pe_end)
723  deallocate(layout2d, global_indices)
724  deallocate(tile1, tile2, tile_id)
725  deallocate(istart1, iend1, jstart1, jend1)
726  deallocate(istart2, iend2, jstart2, jend2)
727 
728 end subroutine setup_domain
729 
730 ! --------------------------------------------------------------------------------------------------
731 
732 subroutine write_geom(self)
733 
734  type(fv3jedi_geom), intent(in) :: self
735 
736  type(fckit_mpi_comm) :: f_comm
737  character(len=255) :: filename
738  integer :: ncid, xf_dimid, yf_dimid, xv_dimid, yv_dimid, ti_dimid, pe_dimid
739  integer :: mydims(3,3), ijdims(1), ijdimf(1), tmpij(1)
740  integer :: varid(8)
741 
742 
743  ! Pointer to fv3jedi geom communicator
744  f_comm = self%f_comm
745 
746  write(filename,"(A9,I0.4,A4)") 'fv3grid_c', self%npx-1, '.nc4'
747 
748  ! Create and open the file for parallel write
749  call nccheck( nf90_create( trim(filename), ior(nf90_netcdf4, nf90_mpiio), ncid, &
750  comm = f_comm%communicator(), info = mpi_info_null), "nf90_create" )
751 
752  !Dimensions
753  call nccheck ( nf90_def_dim(ncid, 'fxdim', self%npx-1 , xf_dimid), "nf90_def_dim fxdim" )
754  call nccheck ( nf90_def_dim(ncid, 'fydim', self%npy-1 , yf_dimid), "nf90_def_dim fydim" )
755  call nccheck ( nf90_def_dim(ncid, 'vxdim', self%npx , xv_dimid), "nf90_def_dim vxdim" )
756  call nccheck ( nf90_def_dim(ncid, 'vydim', self%npy , yv_dimid), "nf90_def_dim vydim" )
757  call nccheck ( nf90_def_dim(ncid, 'ntile', 6 , ti_dimid), "nf90_def_dim ntile" )
758  call nccheck ( nf90_def_dim(ncid, 'nproc', f_comm%size(), pe_dimid), "nf90_def_dim ntile" )
759 
760  !Define variables
761  call nccheck( nf90_def_var(ncid, "flons", nf90_double, (/ xf_dimid, yf_dimid, ti_dimid /), varid(1)), "nf90_def_var flons" )
762  call nccheck( nf90_put_att(ncid, varid(1), "long_name", "longitude of faces") )
763  call nccheck( nf90_put_att(ncid, varid(1), "units", "degrees_east") )
764 
765  call nccheck( nf90_def_var(ncid, "flats", nf90_double, (/ xf_dimid, yf_dimid, ti_dimid /), varid(2)), "nf90_def_var flats" )
766  call nccheck( nf90_put_att(ncid, varid(2), "long_name", "latitude of faces") )
767  call nccheck( nf90_put_att(ncid, varid(2), "units", "degrees_north") )
768 
769  call nccheck( nf90_def_var(ncid, "vlons", nf90_double, (/ xv_dimid, yv_dimid, ti_dimid /), varid(3)), "nf90_def_var vlons" )
770  call nccheck( nf90_put_att(ncid, varid(3), "long_name", "longitude of vertices") )
771  call nccheck( nf90_put_att(ncid, varid(3), "units", "degrees_east") )
772 
773  call nccheck( nf90_def_var(ncid, "vlats", nf90_double, (/ xv_dimid, yv_dimid, ti_dimid /), varid(4)), "nf90_def_var vlats" )
774  call nccheck( nf90_put_att(ncid, varid(4), "long_name", "latitude of vertices") )
775  call nccheck( nf90_put_att(ncid, varid(4), "units", "degrees_north") )
776 
777  call nccheck( nf90_def_var(ncid, "isc", nf90_int, (/ pe_dimid /), varid(5)), "nf90_def_var isc" )
778  call nccheck( nf90_put_att(ncid, varid(5), "long_name", "starting index i direction") )
779  call nccheck( nf90_put_att(ncid, varid(5), "units", "1") )
780 
781  call nccheck( nf90_def_var(ncid, "iec", nf90_int, (/ pe_dimid /), varid(6)), "nf90_def_var iec" )
782  call nccheck( nf90_put_att(ncid, varid(6), "long_name", "ending index i direction") )
783  call nccheck( nf90_put_att(ncid, varid(6), "units", "1") )
784 
785  call nccheck( nf90_def_var(ncid, "jsc", nf90_int, (/ pe_dimid /), varid(7)), "nf90_def_var jsc" )
786  call nccheck( nf90_put_att(ncid, varid(7), "long_name", "starting index j direction") )
787  call nccheck( nf90_put_att(ncid, varid(7), "units", "1") )
788 
789  call nccheck( nf90_def_var(ncid, "jec", nf90_int, (/ pe_dimid /), varid(8)), "nf90_def_var jec" )
790  call nccheck( nf90_put_att(ncid, varid(8), "long_name", "ending index j direction") )
791  call nccheck( nf90_put_att(ncid, varid(8), "units", "1") )
792 
793  ! End define mode
794  call nccheck( nf90_enddef(ncid), "nf90_enddef" )
795 
796  ! Write variables
797  mydims(1,1) = 1; mydims(2,1) = self%npx-1
798  mydims(1,2) = 1; mydims(2,2) = self%npy-1
799  mydims(1,3) = self%ntile; mydims(2,3) = 1
800 
801  call nccheck( nf90_put_var( ncid, varid(1), self%grid_lon(self%isc:self%iec,self%jsc:self%jec), &
802  start = mydims(1,:), count = mydims(2,:) ), "nf90_put_var flons" )
803 
804  call nccheck( nf90_put_var( ncid, varid(2), self%grid_lat(self%isc:self%iec,self%jsc:self%jec), &
805  start = mydims(1,:), count = mydims(2,:) ), "nf90_put_var flats" )
806 
807  mydims(1,1) = 1; mydims(2,1) = self%npx
808  mydims(1,2) = 1; mydims(2,2) = self%npy
809  mydims(1,3) = self%ntile; mydims(2,3) = 1
810 
811  call nccheck( nf90_put_var( ncid, varid(3), self%egrid_lon(self%isc:self%iec+1,self%jsc:self%jec+1), &
812  start = mydims(1,:), count = mydims(2,:) ), "nf90_put_var vlons" )
813 
814  call nccheck( nf90_put_var( ncid, varid(4), self%egrid_lat(self%isc:self%iec+1,self%jsc:self%jec+1), &
815  start = mydims(1,:), count = mydims(2,:) ), "nf90_put_var vlats" )
816 
817  ijdims(1) = f_comm%rank()+1
818  ijdimf(1) = 1
819 
820  tmpij = self%isc
821  call nccheck( nf90_put_var( ncid, varid(5), tmpij, start = ijdims, count = ijdimf ), "nf90_put_var isc" )
822 
823  tmpij = self%iec
824  call nccheck( nf90_put_var( ncid, varid(6), tmpij, start = ijdims, count = ijdimf ), "nf90_put_var iec" )
825 
826  tmpij = self%jsc
827  call nccheck( nf90_put_var( ncid, varid(7), tmpij, start = ijdims, count = ijdimf ), "nf90_put_var jsc" )
828 
829  tmpij = self%jec
830  call nccheck( nf90_put_var( ncid, varid(8), tmpij, start = ijdims, count = ijdimf ), "nf90_put_var jec" )
831 
832  ! Close the file
833  call nccheck ( nf90_close(ncid), "nf90_close" )
834 
835 end subroutine write_geom
836 
837 !--------------------------------------------------------------------------------------------------
838 
839 subroutine getverticalcoordlogp(self, vc, npz, psurf)
840  ! returns log(pressure) at mid level of the vertical column with surface prsssure of psurf
841  ! coded using an example from Jeff Whitaker used in GSI ENKF pacakge
842 
843  type(fv3jedi_geom), intent(in) :: self
844  integer, intent(in) :: npz
845  real(kind=kind_real), intent(in) :: psurf
846  real(kind=kind_real), intent(out) :: vc(npz)
847 
848  real :: plevli(npz+1), plevlm
849  real :: rd, cp, kap, kapr, kap1
850  integer :: k
851 
852  ! constants
853  rd = 2.8705e+2
854  cp = 1.0046e+3
855  kap = rd/cp
856  kapr = cp/rd
857  kap1 = kap + 1.0
858 
859  ! compute interface pressure
860  do k=1,self%npz+1
861  plevli(k) = self%ak(k) + self%bk(k)*psurf
862  enddo
863 
864  ! compute presure at mid level and convert it to logp
865  do k=1,self%npz
866  ! phillips vertical interpolation from guess_grids.F90 in GSI
867  ! (used for global model)
868  plevlm = ((plevli(k)**kap1-plevli(k+1)**kap1)/(kap1*(plevli(k)-plevli(k+1))))**kapr
869  vc(k) = - log(plevlm)
870  enddo
871 
872 end subroutine getverticalcoordlogp
873 ! --------------------------------------------------------------------------------------------------
874 
875 end module fv3jedi_geom_mod
fv3jedi_geom_mod::write_geom
subroutine write_geom(self)
Definition: fv3jedi_geom_mod.f90:733
fv3jedi_netcdf_utils_mod
Definition: fv3jedi_netcdf_utils_mod.F90:6
atm
Definition: atm.F90:1
fv3jedi_geom_mod::set_atlas_lonlat
subroutine set_atlas_lonlat(self, afieldset)
Definition: fv3jedi_geom_mod.f90:536
fv3jedi_geom_mod::create
subroutine create(self, conf, comm, fields)
Definition: fv3jedi_geom_mod.f90:121
fv3jedi_constants_mod::rad2deg
real(kind=kind_real), parameter, public rad2deg
Definition: fv3jedi_constants_mod.f90:13
fv3jedi_geom_mod::delete
subroutine delete(self)
Definition: fv3jedi_geom_mod.f90:478
fv_init_mod
Definition: fv_init_control_mod.f90:6
fv3jedi_geom_mod
Fortran module handling geometry for the FV3 model.
Definition: fv3jedi_geom_mod.f90:8
fv3jedi_constants_mod::ps
real(kind=kind_real), parameter, public ps
Definition: fv3jedi_constants_mod.f90:44
fields_metadata_mod::field_metadata
Definition: fields_metadata_mod.f90:28
fields_metadata_mod
Definition: fields_metadata_mod.f90:6
fv3jedi_geom_mod::getverticalcoordlogp
subroutine, public getverticalcoordlogp(self, vc, npz, psurf)
Definition: fv3jedi_geom_mod.f90:840
fv3jedi_geom_mod::setup_domain
subroutine setup_domain(domain, nx, ny, ntiles, layout_in, io_layout, halo)
Definition: fv3jedi_geom_mod.f90:591
fv3jedi_geom_mod::fv3jedi_geom
Fortran derived type to hold geometry data for the FV3JEDI model.
Definition: fv3jedi_geom_mod.f90:46
fv3jedi_netcdf_utils_mod::nccheck
subroutine, public nccheck(status, iam)
Definition: fv3jedi_netcdf_utils_mod.F90:21
fv3jedi_geom_mod::initialize
subroutine, public initialize(conf, comm)
Definition: fv3jedi_geom_mod.f90:100
fv3jedi_geom_mod::clone
subroutine clone(self, other, fields)
Definition: fv3jedi_geom_mod.f90:356
fv3jedi_constants_mod
Definition: fv3jedi_constants_mod.f90:6
fields_metadata_mod::fields_metadata
Definition: fields_metadata_mod.f90:21
fv3jedi_kinds_mod::kind_real
integer, parameter, public kind_real
Definition: fv3jedi_kinds_mod.f90:14
fv3jedi_geom_mod::fill_atlas_fieldset
subroutine fill_atlas_fieldset(self, afieldset)
Definition: fv3jedi_geom_mod.f90:557
fv3jedi_kinds_mod
Definition: fv3jedi_kinds_mod.f90:6
fv_init_mod::fv_init
subroutine, public fv_init(Atm, dt_atmos_in, grids_on_this_pe, p_split, gtile)
Definition: fv_init_control_mod.f90:28