Line data Source code
1 : module dyn_grid
2 : !-----------------------------------------------------------------------
3 : !
4 : ! Define dynamics computational grid and decomposition.
5 : !
6 : ! Original code: John Drake and Patrick Worley
7 : !
8 : !-----------------------------------------------------------------------
9 :
10 : use shr_kind_mod, only: r8=>shr_kind_r8
11 : use pmgrid, only: plon, plat, plev, plevp, splon, spmd_on
12 : use commap, only: w, w_staggered, clat, clat_staggered, clon, &
13 : latdeg, londeg, latdeg_st, londeg_st
14 : use constituents, only: pcnst
15 : use physconst, only: pi, rearth, omega, spval
16 : use spmd_utils, only: iam
17 : use spmd_dyn, only: spmdinit_dyn, proc, lonrangexy, latrangexy
18 : use time_manager, only: get_step_size
19 :
20 : use pio, only: file_desc_t
21 : use cam_initfiles, only: initial_file_get_id
22 :
23 : use dynamics_vars, only: t_fvdycore_state, t_fvdycore_grid, grid_vars_init
24 : use dyn_internal_state, only: get_dyn_state
25 :
26 : use cam_abortutils, only: endrun
27 : use cam_logfile, only: iulog
28 :
29 : implicit none
30 : private
31 : save
32 :
33 : public :: &
34 : dyn_grid_init, &
35 : dyn_grid_find_gcols, &! find nearest column for given lat/lon
36 : dyn_grid_get_colndx, &! global lat and lon coordinate and MPI process indices
37 : ! corresponding to a specified global column index
38 : dyn_grid_get_elem_coords, &! coordinates of a specified element (latitude)
39 : ! of the dynamics grid (lat slice of the block)
40 : get_block_bounds_d, &! first and last indices in global block ordering
41 : get_block_gcol_d, &! global column indices for given block
42 : get_block_gcol_cnt_d, &! number of columns in given block
43 : get_block_levels_d, &! vertical levels in column
44 : get_block_lvl_cnt_d, &! number of vertical levels in column
45 : get_block_owner_d, &! process "owning" given block
46 : get_dyn_grid_parm, &
47 : get_dyn_grid_parm_real1d, &
48 : get_gcol_block_d, &! global block indices and local columns
49 : ! index for given global column index
50 : get_gcol_block_cnt_d, &! number of blocks containing data
51 : ! from a given global column index
52 : get_horiz_grid_d, &! horizontal grid coordinates
53 : get_horiz_grid_dim_d, &! horizontal dimensions of dynamics grid
54 : physgrid_copy_attributes_d
55 :
56 : ! The FV dynamics grids
57 : integer, parameter, public :: dyn_decomp = 101
58 : integer, parameter, public :: dyn_stagger_decomp = 102 !Backward compatibility
59 : integer, parameter, public :: dyn_ustag_decomp = 102
60 : integer, parameter, public :: dyn_vstag_decomp = 103
61 : integer, parameter, public :: dyn_zonal_decomp = 104
62 :
63 : integer, parameter, public :: ptimelevels = 2 ! number of time levels in the dycore
64 :
65 : integer :: ngcols_d = 0 ! number of dynamics columns
66 :
67 : type(t_fvdycore_grid), pointer :: grid
68 :
69 : !========================================================================================
70 : contains
71 : !========================================================================================
72 :
73 3072 : subroutine dyn_grid_init()
74 :
75 : ! Initialize FV grid, decomposition, and PILGRIM communications
76 :
77 : use hycoef, only: hycoef_init, hyai, hybi, hypi, hypm, nprlev
78 : use ref_pres, only: ref_pres_init
79 :
80 : ! Local variables
81 : type(t_fvdycore_state), pointer :: state
82 :
83 : type(file_desc_t), pointer :: fh_ini
84 :
85 : integer :: i, k, lat
86 :
87 : real(r8) :: dt
88 : real(r8) :: dp
89 : real(r8) :: sum
90 :
91 : character(len=*), parameter :: sub='dyn_grid_init'
92 : !-----------------------------------------------------------------------
93 :
94 : ! Assign pointer to FV internal state object
95 1536 : state => get_dyn_state()
96 : ! Assign pointer to grid object stored in state object
97 1536 : grid => state%grid
98 :
99 : ! Get file handle for initial file and first consistency check
100 1536 : fh_ini => initial_file_get_id()
101 :
102 : ! Set grid size parameters
103 1536 : grid%im = plon
104 1536 : grid%jm = plat
105 1536 : grid%km = plev
106 1536 : grid%kmax = plev + 1
107 1536 : grid%nq = pcnst
108 1536 : grid%ntotq = pcnst
109 :
110 : ! Initialize hybrid coordinate arrays
111 1536 : call hycoef_init(fh_ini)
112 :
113 : ! Initialize reference pressures
114 1536 : call ref_pres_init(hypi, hypm, nprlev)
115 :
116 : ! Hybrid coordinate info for FV grid object
117 1536 : allocate(grid%ak(plev+1), grid%bk(plev+1))
118 1536 : grid%ks = plev
119 110592 : do k = 1, plev+1
120 109056 : grid%ak(k) = hyai(k) * 1.e5_r8
121 109056 : grid%bk(k) = hybi(k)
122 110592 : if (grid%bk(k) == 0._r8) grid%ks = k-1
123 : end do
124 1536 : grid%ptop = grid%ak(1)
125 1536 : grid%pint = grid%ak(grid%ks+1)
126 :
127 : ! Initialize the grid decomposition and PILGRIM communications
128 1536 : call spmdinit_dyn(state%jord, grid)
129 :
130 : ! Initialize FV specific grid object variables
131 1536 : dt = get_step_size()
132 : call grid_vars_init(pi, rearth, omega, dt, state%fft_flt, &
133 1536 : state%am_geom_crrct, grid)
134 :
135 : ! initialize commap variables
136 :
137 : ! latitudes for cell centers
138 1536 : dp = 180._r8/(plat-1)
139 296448 : do lat = 1, plat
140 294912 : latdeg(lat) = -90._r8 + (lat-1)*dp
141 296448 : clat(lat) = latdeg(lat)*pi/180._r8
142 : end do
143 :
144 : ! latitudes for the staggered grid
145 294912 : do lat = 1, plat-1
146 293376 : clat_staggered(lat) = (clat(lat) + clat(lat+1)) / 2._r8
147 294912 : latdeg_st (lat) = clat_staggered(lat)*180._r8/pi
148 : end do
149 :
150 : ! Weights are defined as cos(phi)*(delta-phi)
151 : ! For a sanity check, the sum of w across all lats should be 2.
152 293376 : do lat = 2, plat-1
153 293376 : w(lat) = sin(clat_staggered(lat)) - sin(clat_staggered(lat-1))
154 : end do
155 1536 : w(1) = sin(clat_staggered(1)) + 1._r8
156 1536 : w(plat) = w(1)
157 :
158 1536 : sum = 0._r8
159 296448 : do lat=1,plat
160 296448 : sum = sum + w(lat)
161 : end do
162 1536 : if (abs(sum - 2._r8) > 1.e-8_r8) then
163 0 : write(iulog,*) sub//': ERROR: weights do not sum to 2. sum=', sum
164 0 : call endrun(sub//': ERROR: weights do not sum to 2.')
165 : end if
166 :
167 294912 : dp = pi / real(plat-1,r8)
168 294912 : do lat = 1, plat-1
169 294912 : w_staggered(lat) = sin(clat(lat+1)) - sin(clat(lat))
170 : end do
171 :
172 1536 : sum = 0._r8
173 294912 : do lat = 1, plat-1
174 294912 : sum = sum + w_staggered(lat)
175 : end do
176 :
177 1536 : if (abs(sum - 2._r8) > 1.e-8_r8) then
178 0 : write(iulog,*) sub//': ERROR: staggered weights do not sum to 2. sum=', sum
179 0 : call endrun(sub//': ERROR: staggered weights do not sum to 2.')
180 : end if
181 :
182 : ! longitudes for cell centers
183 296448 : do lat = 1, plat
184 85231104 : do i = 1, plon
185 84934656 : londeg(i,lat) = (i-1)*360._r8/plon
186 85229568 : clon(i,lat) = (i-1)*2._r8*pi/plon
187 : end do
188 : end do
189 :
190 : ! longitudes for staggered grid
191 296448 : do lat = 1, plat
192 85231104 : do i = 1, splon
193 85229568 : londeg_st(i,lat) = (i-1.5_r8)*360._r8/splon
194 : end do
195 : end do
196 :
197 : ! Define the CAM grids
198 1536 : call define_cam_grids()
199 :
200 1536 : end subroutine dyn_grid_init
201 :
202 : !========================================================================================
203 :
204 3072 : subroutine get_block_bounds_d(block_first, block_last)
205 :
206 : ! Return first and last indices used in global block ordering
207 :
208 : ! Arguments
209 : integer, intent(out) :: block_first ! first (global) index used for blocks
210 : integer, intent(out) :: block_last ! last (global) index used for blocks
211 : !---------------------------------------------------------------------------
212 :
213 3072 : block_first = 1
214 3072 : if (spmd_on .eq. 1) then
215 : ! Assume 1 block per subdomain
216 3072 : block_last = grid%nprxy_x*grid%nprxy_y
217 : else
218 : ! latitude slice block
219 0 : block_last = plat
220 : end if
221 :
222 1536 : end subroutine get_block_bounds_d
223 :
224 : !========================================================================================
225 :
226 1179648 : subroutine get_block_gcol_d(blockid, size, cdex)
227 :
228 : ! Return list of dynamics column indices in given block
229 :
230 : ! Arguments
231 : integer, intent(in) :: blockid ! global block id
232 : integer, intent(in) :: size ! array size
233 :
234 : integer, intent(out):: cdex(size) ! global column indices
235 :
236 : ! Local workspace
237 : integer :: i,j ! block coordinates
238 : integer :: blksiz ! block size
239 : integer :: k,l ! loop indices
240 : integer :: n ! column index
241 : character(len=*), parameter :: sub='get_block_gcol_d'
242 : !---------------------------------------------------------------------------
243 :
244 1179648 : if (spmd_on .eq. 1) then
245 1179648 : j = (blockid-1) / grid%nprxy_x + 1
246 1179648 : i = blockid - (j-1) * grid%nprxy_x
247 : #if ( defined SPMD )
248 1179648 : blksiz = (lonrangexy(2,i)-lonrangexy(1,i)+1) * &
249 2359296 : (latrangexy(2,j)-latrangexy(1,j)+1)
250 1179648 : if (size < blksiz) then
251 0 : write(iulog,*) sub//': ERROR: array not large enough (', &
252 0 : size,' < ',blksiz,' ) '
253 0 : call endrun(sub//': ERROR: array not large enough')
254 : else
255 1179648 : n = 0
256 4718592 : do k=latrangexy(1,j),latrangexy(2,j)
257 89653248 : do l=lonrangexy(1,i),lonrangexy(2,i)
258 84934656 : n = n + 1
259 88473600 : cdex(n) = l + (k-1)*plon
260 : end do
261 : end do
262 : end if
263 : #endif
264 : else
265 0 : if (size < plon) then
266 0 : write(iulog,*)'GET_BLOCK_GCOL_D: array not large enough (', &
267 0 : size,' < ',plon,' ) '
268 0 : call endrun
269 : else
270 0 : n = (blockid-1)*plon
271 0 : do i = 1,plon
272 0 : n = n + 1
273 0 : cdex(i) = n
274 : end do
275 : end if
276 : end if
277 :
278 1179648 : end subroutine get_block_gcol_d
279 :
280 : !========================================================================================
281 :
282 2360832 : integer function get_block_gcol_cnt_d(blockid)
283 :
284 : ! Return number of dynamics columns in indicated block
285 :
286 : ! Arguments
287 : integer, intent(in) :: blockid ! global block id
288 :
289 : ! Local workspace
290 : integer :: i, j
291 : !---------------------------------------------------------------------------
292 :
293 2360832 : if (spmd_on .eq. 1) then
294 2360832 : j = (blockid-1) / grid%nprxy_x + 1
295 2360832 : i = blockid - (j-1) * grid%nprxy_x
296 : #if ( defined SPMD )
297 2360832 : get_block_gcol_cnt_d = (lonrangexy(2,i)-lonrangexy(1,i)+1) * &
298 4721664 : (latrangexy(2,j)-latrangexy(1,j)+1)
299 : #endif
300 : else
301 : get_block_gcol_cnt_d = plon
302 : end if
303 :
304 2360832 : end function get_block_gcol_cnt_d
305 :
306 : !========================================================================================
307 :
308 112128 : integer function get_block_lvl_cnt_d(blockid, bcid)
309 :
310 : ! Return number of levels in indicated column. If column
311 : ! includes surface fields, then it is defined to also
312 : ! include level 0.
313 :
314 : ! Arguments
315 : integer, intent(in) :: blockid ! global block id
316 : integer, intent(in) :: bcid ! column index within block
317 : !-----------------------------------------------------------------------
318 :
319 : ! latitude slice block
320 112128 : get_block_lvl_cnt_d = plev + 1
321 :
322 112128 : end function get_block_lvl_cnt_d
323 :
324 : !========================================================================================
325 :
326 110592 : subroutine get_block_levels_d(blockid, bcid, lvlsiz, levels)
327 :
328 : ! Return level indices in indicated column. If column
329 : ! includes surface fields, then it is defined to also
330 : ! include level 0.
331 :
332 : ! Arguments
333 : integer, intent(in) :: blockid ! global block id
334 : integer, intent(in) :: bcid ! column index within block
335 : integer, intent(in) :: lvlsiz ! dimension of levels array
336 :
337 : integer, intent(out) :: levels(lvlsiz) ! levels indices for block
338 :
339 : ! Local workspace
340 : integer :: k ! loop index
341 : character(len=*), parameter :: sub='get_block_levels_d'
342 : !-----------------------------------------------------------------------
343 :
344 : ! latitude slice block
345 110592 : if (lvlsiz < plev + 1) then
346 0 : write(iulog,*) sub//': ERROR: levels array not large enough (', &
347 0 : lvlsiz,' < ',plev + 1,' ) '
348 0 : call endrun(sub//': ERROR: levels array not large enough')
349 : else
350 7962624 : do k = 0, plev
351 7962624 : levels(k+1) = k
352 : end do
353 110592 : do k = plev+2, lvlsiz
354 110592 : levels(k) = -1
355 : end do
356 : end if
357 :
358 110592 : end subroutine get_block_levels_d
359 :
360 : !========================================================================================
361 :
362 465960960 : subroutine get_gcol_block_d(gcol, cnt, blockid, bcid, localblockid)
363 :
364 : ! Return global block index and local column index
365 : ! for global column index
366 :
367 : ! Arguments
368 : integer, intent(in) :: gcol ! global column index
369 : integer, intent(in) :: cnt ! size of blockid and bcid arrays
370 :
371 : integer, intent(out) :: blockid(cnt) ! block index
372 : integer, intent(out) :: bcid(cnt) ! column index within block
373 : integer, intent(out), optional :: localblockid(cnt)
374 :
375 : ! Local workspace
376 : integer :: i, j, ii, jj ! loop indices
377 : integer :: glon, glat ! global longitude and latitude indices
378 : integer :: ddlon ! number of longitudes in block
379 : character(len=*), parameter :: sub='get_gcol_block_d'
380 : !---------------------------------------------------------------------------
381 :
382 : ! lon/lat block
383 465960960 : if (cnt < 1) then
384 0 : write(iulog,*) sub//': ERROR: arrays not large enough (', cnt,' < ',1,' )'
385 0 : call endrun(sub//': ERROR: arrays not large enough')
386 : else
387 465960960 : if (spmd_on == 1) then
388 : ! Determine global latitude and longitude coordinate indices from
389 : ! global column index
390 465960960 : glat = (gcol-1)/plon + 1
391 465960960 : glon = gcol - ((glat-1)*plon)
392 :
393 : ! Determine block coordinates (ii,jj), where ii ranges from 1 to
394 : ! nprxy_x and jj ranges from 1 to nprxy_y.
395 : #if ( defined SPMD )
396 465960960 : ii = 0
397 6057492480 : do i = 1, grid%nprxy_x
398 6057492480 : if (lonrangexy(1,i) <= glon .and. glon <= lonrangexy(2,i)) ii=i
399 : end do
400 465960960 : jj = 0
401 30287462400 : do j = 1, grid%nprxy_y
402 30287462400 : if (latrangexy(1,j) <= glat .and. glat <= latrangexy(2,j)) jj=j
403 : end do
404 465960960 : if (ii == 0 .or. jj == 0) then
405 0 : write(iulog,*) sub//': ERROR: could not find block indices for (', &
406 0 : glon,',',glat,' ) '
407 0 : call endrun(sub//': ERROR: could not find block indices')
408 : end if
409 :
410 : ! Global block index
411 465960960 : blockid(1) = (jj-1)*grid%nprxy_x+ii
412 :
413 : ! Local coordinates in block
414 465960960 : j = glat-latrangexy(1,jj)+1
415 465960960 : i = glon-lonrangexy(1,ii)+1
416 465960960 : ddlon = lonrangexy(2,ii)-lonrangexy(1,ii)+1
417 :
418 : ! Local column index in block
419 465960960 : bcid(1) = (j-1)*ddlon+i
420 : #endif
421 : else
422 0 : glat = (gcol-1)/plon + 1
423 0 : glon = gcol - ((glat-1)*plon)
424 :
425 0 : blockid(1) = glat
426 0 : bcid(1) = glon
427 : end if
428 :
429 32133611520 : do j=2,cnt
430 31667650560 : blockid(j) = -1
431 32133611520 : bcid(j) = -1
432 : end do
433 :
434 : end if
435 :
436 465960960 : end subroutine get_gcol_block_d
437 :
438 : !========================================================================================
439 :
440 424673280 : integer function get_gcol_block_cnt_d(gcol)
441 :
442 : ! Return number of blocks contain data for the vertical column
443 : ! with the given global column index
444 :
445 : ! Arguments
446 : integer, intent(in) :: gcol ! global column index
447 : !---------------------------------------------------------------------------
448 :
449 : ! lon/lat block
450 424673280 : get_gcol_block_cnt_d = 1
451 :
452 424673280 : end function get_gcol_block_cnt_d
453 :
454 : !========================================================================================
455 :
456 468320256 : integer function get_block_owner_d(blockid)
457 :
458 : ! Return id of processor that "owns" the indicated block
459 :
460 : ! Arguments
461 : integer, intent(in) :: blockid ! global block id
462 : !---------------------------------------------------------------------------
463 :
464 : ! latitude slice block
465 : #if (defined SPMD)
466 468320256 : if (spmd_on .eq. 1) then
467 468320256 : get_block_owner_d = blockid - 1
468 : else
469 0 : get_block_owner_d = proc(blockid)
470 : end if
471 : #else
472 : get_block_owner_d = 0
473 : #endif
474 :
475 468320256 : end function get_block_owner_d
476 :
477 : !========================================================================================
478 :
479 6144 : subroutine get_horiz_grid_dim_d(hdim1_d, hdim2_d)
480 :
481 : ! Returns declared horizontal dimensions of computational grid.
482 : ! Note that global column ordering is assumed to be compatible
483 : ! with the first dimension major ordering of the 2D array.
484 :
485 : ! Arguments
486 : integer, intent(out) :: hdim1_d ! first horizontal dimension
487 : integer, intent(out) :: hdim2_d ! second horizontal dimension
488 : !---------------------------------------------------------------------------
489 :
490 6144 : if (ngcols_d == 0) then
491 1536 : ngcols_d = plat*plon
492 : end if
493 6144 : hdim1_d = plon
494 6144 : hdim2_d = plat
495 :
496 6144 : end subroutine get_horiz_grid_dim_d
497 :
498 : !========================================================================================
499 :
500 6144 : subroutine get_horiz_grid_d(size, clat_d_out, clon_d_out, area_d_out, wght_d_out, &
501 3072 : lat_d_out, lon_d_out)
502 :
503 : ! Return latitude and longitude (in radians), column surface
504 : ! area (in radians squared) and surface integration weights
505 : ! for global column indices that will be passed to/from physics
506 :
507 : ! Arguments
508 : integer, intent(in) :: size ! array sizes
509 :
510 : real(r8), intent(out), optional :: clat_d_out(size) ! column latitudes
511 : real(r8), intent(out), optional :: clon_d_out(size) ! column longitudes
512 :
513 : real(r8), intent(out), optional :: area_d_out(size) ! column surface
514 : ! area
515 : real(r8), intent(out), optional :: wght_d_out(size) ! column integration
516 : ! weight
517 : real(r8), intent(out), optional :: lat_d_out(size) ! column deg latitudes
518 : real(r8), intent(out), optional :: lon_d_out(size) ! column deg longitudes
519 :
520 : ! Local workspace
521 : integer :: i, j ! loop indices
522 : integer :: n ! column index
523 : real(r8) :: ns_vert(2,plon) ! latitude grid vertices
524 : real(r8) :: ew_vert(2,plon) ! longitude grid vertices
525 : real(r8) :: del_theta ! difference in latitude at a grid point
526 : real(r8) :: del_phi ! difference in longitude at a grid point
527 : real(r8), parameter :: degtorad=pi/180.0_r8 ! convert degrees to radians
528 : character(len=128) :: errormsg
529 : character(len=*), parameter :: sub='get_horiz_grid_d'
530 : !----------------------------------------------------------------------------
531 :
532 3072 : if (present(clon_d_out)) then
533 1536 : if (size == ngcols_d) then
534 : n = 0
535 296448 : do j = 1, plat
536 85231104 : do i = 1, plon
537 84934656 : n = n + 1
538 85229568 : clon_d_out(n) = clon(i,j)
539 : end do
540 : end do
541 0 : else if(size == plon) then
542 0 : clon_d_out(:) = clon(:,1)
543 : else
544 0 : write(errormsg, '(a,4(i0,a))')'clon_d_out array size incorrect (', &
545 0 : size, ' /= ', ngcols_d, ' .and. ', size, ' /= ', plon,') '
546 0 : call endrun(sub//': ERROR: '//errormsg)
547 : end if
548 : end if
549 :
550 3072 : if (present(clat_d_out)) then
551 1536 : if (size == ngcols_d) then
552 : n = 0
553 296448 : do j = 1, plat
554 85231104 : do i = 1, plon
555 84934656 : n = n + 1
556 85229568 : clat_d_out(n) = clat(j)
557 : end do
558 : end do
559 0 : else if (size == plat) then
560 0 : clat_d_out(:) = clat(:)
561 : else
562 0 : write(errormsg, '(a,4(i0,a))')'clat_d_out array size incorrect (', &
563 0 : size, ' /= ', ngcols_d, ' .and. ', size, ' /= ', plat,') '
564 0 : call endrun(sub//': ERROR: '//errormsg)
565 : end if
566 : end if
567 :
568 3072 : if (size==plat .and. present(wght_d_out)) then
569 0 : wght_d_out(:) = w(:)
570 0 : return
571 : end if
572 :
573 3072 : if ( ( present(area_d_out) ) .or. ( present(wght_d_out) ) ) then
574 1536 : if ((size < ngcols_d) .and. present(area_d_out)) then
575 0 : write(errormsg, '(a,2(i0,a))')'area_d_out array size incorrect (', &
576 0 : size, ' /= ', ngcols_d, ') '
577 0 : call endrun(sub//': ERROR: '//errormsg)
578 : else if ((size < ngcols_d) .and. present(area_d_out)) then
579 : write(errormsg, '(a,2(i0,a))')'wght_d_out array size incorrect (', &
580 : size, ' /= ', ngcols_d, ') '
581 : call endrun(sub//': ERROR: '//errormsg)
582 : end if
583 :
584 : n = 0
585 296448 : do j = 1, plat
586 :
587 : ! First, determine vertices of each grid point.
588 : ! Verticies are ordered as follows:
589 : ! ns_vert: 1=lower left, 2 = upper left
590 : ! ew_vert: 1=lower left, 2 = lower right
591 :
592 : ! Latitude vertices
593 255098880 : ns_vert(:,:) = spval
594 294912 : if (j .eq. 1) then
595 443904 : ns_vert(1,:plon) = -90._r8 + (latdeg(1) - latdeg(2))*0.5_r8
596 : else
597 84785664 : ns_vert(1,:plon) = (latdeg(j) + latdeg(j-1) )*0.5_r8
598 : end if
599 :
600 294912 : if (j .eq. plat) then
601 443904 : ns_vert(2,:plon) = 90._r8 + (latdeg(plat) - latdeg(plat-1))*0.5_r8
602 : else
603 84785664 : ns_vert(2,:plon) = (latdeg(j) + latdeg(j+1) )*0.5_r8
604 : end if
605 :
606 : ! Longitude vertices
607 255098880 : ew_vert(:,:) = spval
608 294912 : ew_vert(1,1) = (londeg(1,j) - 360._r8 + londeg(plon,j))*0.5_r8
609 84934656 : ew_vert(1,2:plon) = (londeg(1:plon-1,j)+ londeg(2:plon,j))*0.5_r8
610 84934656 : ew_vert(2,:plon-1) = ew_vert(1,2:plon)
611 294912 : ew_vert(2,plon) = (londeg(plon,j) + (360._r8 + londeg(1,j)))*0.5_r8
612 :
613 85231104 : do i = 1,plon
614 84934656 : n = n + 1
615 :
616 84934656 : if (j .eq. 1) then
617 442368 : del_phi = -sin( latdeg(j)*degtorad ) + sin( ns_vert(2,i)*degtorad )
618 84492288 : else if (j .eq. plat) then
619 442368 : del_phi = sin( latdeg(j)*degtorad ) - sin( ns_vert(1,i)*degtorad )
620 : else
621 84049920 : del_phi = sin( ns_vert(2,i)*degtorad ) - sin( ns_vert(1,i)*degtorad )
622 : end if
623 :
624 84934656 : del_theta = ( ew_vert(2,i) - ew_vert(1,i) )*degtorad
625 :
626 84934656 : if (present(area_d_out)) area_d_out(n) = del_theta*del_phi
627 85229568 : if (present(wght_d_out)) wght_d_out(n) = del_theta*del_phi
628 : end do
629 : end do
630 : end if
631 :
632 3072 : if (present(lon_d_out)) then
633 1536 : if (size == ngcols_d) then
634 : n = 0
635 296448 : do j = 1, plat
636 85231104 : do i = 1, plon
637 84934656 : n = n + 1
638 85229568 : lon_d_out(n) = londeg(i,j)
639 : end do
640 : end do
641 0 : else if(size == plon) then
642 0 : lon_d_out(:) = londeg(:,1)
643 : else
644 0 : write(errormsg, '(a,4(i0,a))')'lon_d_out array size incorrect (', &
645 0 : size, ' /= ', ngcols_d, ' .and. ', size, ' /= ', plon,') '
646 0 : call endrun(sub//': ERROR: '//errormsg)
647 : end if
648 : end if
649 :
650 3072 : if (present(lat_d_out)) then
651 1536 : if (size == ngcols_d) then
652 : n = 0
653 296448 : do j = 1, plat
654 85231104 : do i = 1, plon
655 84934656 : n = n + 1
656 85229568 : lat_d_out(n) = latdeg(j)
657 : end do
658 : end do
659 0 : else if (size == plat) then
660 0 : lat_d_out(:) = latdeg(:)
661 : else
662 0 : write(errormsg, '(a,4(i0,a))')'lat_d_out array size incorrect (', &
663 0 : size, ' /= ', ngcols_d, ' .and. ', size, ' /= ', plat,') '
664 0 : call endrun(sub//': ERROR: '//errormsg)
665 : end if
666 : end if
667 :
668 10752 : end subroutine get_horiz_grid_d
669 :
670 : !========================================================================================
671 :
672 : function get_dyn_grid_parm_real2d(name) result(rval)
673 :
674 : character(len=*), intent(in) :: name
675 : real(r8), pointer :: rval(:,:)
676 :
677 : if(name.eq.'clon') then
678 : rval => clon
679 : else if(name.eq.'londeg') then
680 : rval => londeg
681 : else if(name.eq.'londeg_st') then
682 : rval => londeg_st
683 : else
684 : nullify(rval)
685 : end if
686 : end function get_dyn_grid_parm_real2d
687 :
688 : !========================================================================================
689 :
690 0 : function get_dyn_grid_parm_real1d(name) result(rval)
691 :
692 : character(len=*), intent(in) :: name
693 : real(r8), pointer :: rval(:)
694 :
695 0 : if(name.eq.'clat') then
696 0 : rval => clat
697 0 : else if(name.eq.'latdeg') then
698 0 : rval => latdeg
699 0 : else if(name.eq.'latdeg_st') then
700 0 : rval => latdeg_st
701 0 : else if(name.eq.'clatdeg_staggered') then
702 0 : rval => latdeg_st
703 0 : else if(name.eq.'w') then
704 0 : rval => w
705 0 : else if(name.eq.'w_staggered') then
706 0 : rval => w_staggered
707 : else
708 0 : nullify(rval)
709 : end if
710 0 : end function get_dyn_grid_parm_real1d
711 :
712 : !========================================================================================
713 :
714 270336 : integer function get_dyn_grid_parm(name) result(ival)
715 :
716 : character(len=*), intent(in) :: name
717 :
718 270336 : if (name.eq.'splon') then
719 : ival = splon
720 270336 : else if (name.eq.'beglonxy') then
721 4608 : ival = grid%ifirstxy
722 265728 : else if (name.eq.'endlonxy') then
723 4608 : ival = grid%ilastxy
724 261120 : else if (name.eq.'beglatxy') then
725 4608 : ival = grid%jfirstxy
726 256512 : else if (name.eq.'endlatxy') then
727 4608 : ival = grid%jlastxy
728 251904 : else if (name.eq.'plat') then
729 : ival = plat
730 125952 : else if (name.eq.'plon') then
731 : ival = plon
732 0 : else if (name.eq.'plev') then
733 : ival = plev
734 0 : else if (name.eq.'plevp') then
735 : ival = plevp
736 : else
737 0 : ival = -1
738 : end if
739 :
740 270336 : end function get_dyn_grid_parm
741 :
742 : !========================================================================================
743 :
744 0 : subroutine dyn_grid_find_gcols(lat, lon, nclosest, owners, indx, &
745 0 : jndx, rlat, rlon, idyn_dists)
746 :
747 : ! Return the lat/lon information (and corresponding MPI task numbers (owners))
748 : ! of the global model grid columns nearest to the input coordinate (lat,lon)
749 :
750 : ! arguments
751 : real(r8), intent(in) :: lat
752 : real(r8), intent(in) :: lon
753 : integer, intent(in) :: nclosest
754 : integer, intent(out) :: owners(nclosest)
755 : integer, intent(out) :: indx(nclosest)
756 : integer, intent(out) :: jndx(nclosest)
757 :
758 : real(r8),optional, intent(out) :: rlon(nclosest)
759 : real(r8),optional, intent(out) :: rlat(nclosest)
760 : real(r8),optional, intent(out) :: idyn_dists(nclosest)
761 :
762 : ! local variables
763 : real(r8) :: dist ! the distance (in radians**2 from lat, lon)
764 : real(r8) :: latr, lonr ! lat, lon inputs converted to radians
765 : integer :: ngcols
766 : integer :: i, j
767 :
768 : integer :: blockid(1), bcid(1), lclblockid(1)
769 :
770 0 : real(r8), allocatable :: clat_d(:), clon_d(:), distmin(:)
771 0 : integer, allocatable :: igcol(:)
772 : real(r8), parameter :: rad2deg = 180._r8/pi
773 : !----------------------------------------------------------------------------
774 :
775 0 : latr = lat/rad2deg
776 0 : lonr = lon/rad2deg
777 :
778 0 : ngcols = plon*plat
779 0 : allocate( clat_d(1:ngcols) )
780 0 : allocate( clon_d(1:ngcols) )
781 0 : allocate( igcol(nclosest) )
782 0 : allocate( distmin(nclosest) )
783 :
784 0 : call get_horiz_grid_d(ngcols, clat_d_out=clat_d, clon_d_out=clon_d)
785 :
786 0 : igcol(:) = -999
787 0 : distmin(:) = 1.e10_r8
788 :
789 0 : do i = 1, ngcols
790 :
791 : ! Use the Spherical Law of Cosines to find the great-circle distance.
792 0 : dist = acos(sin(latr) * sin(clat_d(i)) + cos(latr) * cos(clat_d(i)) * &
793 0 : cos(clon_d(i) - lonr)) * rearth
794 0 : do j = nclosest, 1, -1
795 0 : if (dist < distmin(j)) then
796 :
797 0 : if (j < nclosest) then
798 0 : distmin(j+1) = distmin(j)
799 0 : igcol(j+1) = igcol(j)
800 : end if
801 :
802 0 : distmin(j) = dist
803 0 : igcol(j) = i
804 : else
805 : exit
806 : end if
807 : end do
808 : end do
809 :
810 0 : do i = 1,nclosest
811 :
812 0 : call get_gcol_block_d( igcol(i), 1, blockid, bcid, lclblockid )
813 0 : owners(i) = get_block_owner_d(blockid(1))
814 :
815 0 : if ( iam==owners(i) ) then
816 : ! get global lat and lon coordinate indices from global column index
817 : ! -- plon is global number of longitude grid points
818 0 : jndx(i) = (igcol(i)-1)/plon + 1
819 0 : indx(i) = igcol(i) - (jndx(i)-1)*plon
820 : else
821 0 : jndx(i) = -1
822 0 : indx(i) = -1
823 : endif
824 :
825 0 : if ( present(rlat) ) rlat(i) = clat_d(igcol(i)) * rad2deg
826 0 : if ( present(rlon) ) rlon(i) = clon_d(igcol(i)) * rad2deg
827 :
828 0 : if (present(idyn_dists)) then
829 0 : idyn_dists(i) = distmin(i)
830 : end if
831 :
832 : end do
833 :
834 0 : deallocate( clat_d )
835 0 : deallocate( clon_d )
836 0 : deallocate( igcol )
837 0 : deallocate( distmin )
838 :
839 0 : end subroutine dyn_grid_find_gcols
840 :
841 : !========================================================================================
842 :
843 0 : subroutine dyn_grid_get_colndx(igcol, nclosest, owners, indx, jndx)
844 :
845 : ! arguments
846 : integer, intent(in) :: nclosest
847 : integer, intent(in) :: igcol(nclosest)
848 : integer, intent(out) :: owners(nclosest)
849 : integer, intent(out) :: indx(nclosest)
850 : integer, intent(out) :: jndx(nclosest)
851 :
852 : integer :: i
853 : integer :: blockid(1), bcid(1), lclblockid(1)
854 :
855 0 : do i = 1,nclosest
856 :
857 0 : call get_gcol_block_d(igcol(i), 1, blockid, bcid, lclblockid)
858 0 : owners(i) = get_block_owner_d(blockid(1))
859 :
860 0 : if ( iam==owners(i) ) then
861 : ! get global lat and lon coordinate indices from global column index
862 : ! -- plon is global number of longitude grid points
863 0 : jndx(i) = (igcol(i)-1)/plon + 1
864 0 : indx(i) = igcol(i) - (jndx(i)-1)*plon
865 : else
866 0 : jndx(i) = -1
867 0 : indx(i) = -1
868 : end if
869 :
870 : end do
871 :
872 0 : end subroutine dyn_grid_get_colndx
873 :
874 : !========================================================================================
875 :
876 0 : subroutine dyn_grid_get_elem_coords( latndx, rlon, rlat, cdex )
877 :
878 : ! return coordinates of a latitude slice of the block corresponding
879 : ! to latitude index latndx
880 :
881 : ! arguments
882 : integer, intent(in) :: latndx ! lat index
883 :
884 : real(r8),optional, intent(out) :: rlon(:) ! longitudes of the columns in the latndx slice
885 : real(r8),optional, intent(out) :: rlat(:) ! latitudes of the columns in the latndx slice
886 : integer, optional, intent(out) :: cdex(:) ! global column index
887 :
888 : integer :: i, ii, j
889 : !----------------------------------------------------------------------------
890 :
891 0 : if (present(cdex)) cdex(:) = -1
892 0 : if (present(rlat)) rlat(:) = -999._r8
893 0 : if (present(rlon)) rlon(:) = -999._r8
894 :
895 0 : j = latndx
896 0 : ii = 0
897 0 : do i = grid%ifirstxy, grid%ilastxy
898 0 : ii = ii+1
899 0 : if (present(cdex)) cdex(ii) = i + (j-1)*plon
900 0 : if (present(rlat)) rlat(ii) = clat(j)
901 0 : if (present(rlon)) rlon(ii) = clon(i,1)
902 : end do
903 :
904 0 : end subroutine dyn_grid_get_elem_coords
905 :
906 : !========================================================================================
907 :
908 1536 : subroutine define_cam_grids()
909 :
910 : use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap
911 : use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
912 : use cam_grid_support, only: cam_grid_attribute_copy
913 :
914 : integer :: i, j, ind
915 : integer :: beglonxy, endlonxy
916 : integer :: beglatxy, endlatxy
917 : type(horiz_coord_t), pointer :: lat_coord
918 : type(horiz_coord_t), pointer :: lon_coord
919 : type(horiz_coord_t), pointer :: slat_coord
920 : type(horiz_coord_t), pointer :: slon_coord
921 : type(horiz_coord_t), pointer :: zlon_coord
922 1536 : integer(iMap), allocatable :: coord_map(:)
923 1536 : integer(iMap), pointer :: grid_map(:,:)
924 : real(r8) :: zlon_bnds(2,1)
925 1536 : real(r8), allocatable :: latvals(:)
926 1536 : real(r8), pointer :: rattval(:)
927 : logical :: is_lon_distributed
928 : !-----------------------------------------------------------------------------
929 :
930 : ! Note: not using get_horiz_grid_dim_d or get_horiz_grid_d since those
931 : ! are deprecated ('cause I said so' -- goldy)
932 :
933 1536 : nullify(lat_coord)
934 1536 : nullify(lon_coord)
935 1536 : nullify(slat_coord)
936 1536 : nullify(slon_coord)
937 1536 : nullify(zlon_coord)
938 1536 : nullify(grid_map)
939 1536 : nullify(rattval)
940 :
941 1536 : beglonxy = grid%ifirstxy
942 1536 : endlonxy = grid%ilastxy
943 1536 : beglatxy = grid%jfirstxy
944 1536 : endlatxy = grid%jlastxy
945 :
946 1536 : if (iam >= grid%npes_xy) then
947 : ! NB: On inactive PEs, beglonxy should be one and endlonxy should be zero
948 0 : if (beglonxy /= 1) then
949 0 : call endrun("DEFINE_CAM_GRIDS: ERROR: Bad beglonxy")
950 : end if
951 0 : if (endlonxy /= 0) then
952 0 : call endrun("DEFINE_CAM_GRIDS: ERROR: Bad endlonxy")
953 : end if
954 : ! NB: On inactive PEs, beglatxy should be one and endlatxy should be zero
955 0 : if (beglatxy /= 1) then
956 0 : call endrun("DEFINE_CAM_GRIDS: ERROR: Bad beglatxy")
957 : end if
958 0 : if (endlatxy /= 0) then
959 0 : call endrun("DEFINE_CAM_GRIDS: ERROR: Bad endlatxy")
960 : end if
961 : end if
962 :
963 : ! Figure out if lon and slon are distributed dimensions
964 1536 : is_lon_distributed = (grid%nprxy_x > 1)
965 :
966 : ! Grid for cell centers
967 : ! Make a map
968 4608 : allocate(grid_map(4, ((endlonxy - beglonxy + 1) * (endlatxy - beglatxy + 1))))
969 1536 : ind = 0
970 6144 : do i = beglatxy, endlatxy
971 116736 : do j = beglonxy, endlonxy
972 110592 : ind = ind + 1
973 110592 : grid_map(1, ind) = j
974 110592 : grid_map(2, ind) = i
975 110592 : grid_map(3, ind) = j
976 115200 : grid_map(4, ind) = i
977 : end do
978 : end do
979 : ! Cell-centered latitude coordinate
980 4608 : allocate(coord_map(endlatxy - beglatxy + 1))
981 1536 : if (endlatxy >= beglatxy) then
982 1536 : if (beglonxy == 1) then
983 1408 : coord_map = (/ (i, i = beglatxy, endlatxy) /)
984 : else
985 5632 : coord_map = 0
986 : end if
987 : end if
988 : lat_coord => horiz_coord_create('lat', '', plat, 'latitude', &
989 0 : 'degrees_north', beglatxy, endlatxy, latdeg(beglatxy:endlatxy), &
990 1536 : map=coord_map)
991 1536 : deallocate(coord_map)
992 :
993 : ! Cell-centered longitude coordinate
994 1536 : if (is_lon_distributed) then
995 4608 : allocate(coord_map(endlonxy - beglonxy + 1))
996 1536 : if (endlonxy >= beglonxy) then
997 1536 : if (beglatxy == 1) then
998 1776 : coord_map = (/ (i, i = beglonxy, endlonxy) /)
999 : else
1000 37800 : coord_map = 0
1001 : end if
1002 : end if
1003 : lon_coord => horiz_coord_create('lon', '', plon, 'longitude', &
1004 0 : 'degrees_east', beglonxy, endlonxy, londeg(beglonxy:endlonxy,1), &
1005 1536 : map=coord_map)
1006 1536 : deallocate(coord_map)
1007 : else
1008 : lon_coord => horiz_coord_create('lon', '', plon, 'longitude', &
1009 0 : 'degrees_east', beglonxy, endlonxy, londeg(beglonxy:endlonxy,1))
1010 : end if
1011 : ! Cell-centered grid
1012 : call cam_grid_register('fv_centers', dyn_decomp, lat_coord, lon_coord, &
1013 1536 : grid_map, unstruct=.false.)
1014 1536 : allocate(rattval(size(w)))
1015 592896 : rattval = w
1016 1536 : call cam_grid_attribute_register('fv_centers', 'gw', 'latitude weights', 'lat', rattval)
1017 1536 : nullify(rattval)
1018 1536 : nullify(grid_map) ! Belongs to the grid
1019 :
1020 : ! Staggered grid for U_S
1021 : ! Make a map
1022 4608 : allocate(grid_map(4, ((endlonxy - beglonxy + 1) * (endlatxy - beglatxy + 1))))
1023 1536 : ind = 0
1024 6144 : do i = beglatxy, endlatxy
1025 116736 : do j = beglonxy, endlonxy
1026 110592 : ind = ind + 1
1027 110592 : grid_map(1, ind) = j
1028 110592 : grid_map(2, ind) = i
1029 110592 : grid_map(3, ind) = j
1030 115200 : if ((i == beglatxy) .and. (beglatxy == 1)) then
1031 576 : grid_map(4, ind) = 0
1032 : else
1033 110016 : grid_map(4, ind) = i - 1
1034 : end if
1035 : end do
1036 : end do
1037 :
1038 : ! Staggered latitudes 'skip' the first one so they are 'off by one'
1039 : ! This means we always must have a coordinate map
1040 4608 : allocate(coord_map(endlatxy - beglatxy + 1))
1041 : ! NB: coord_map(1) == 0 when beglat == 1, that element is not output
1042 6144 : do i = 1, size(coord_map)
1043 6144 : if (beglonxy == 1) then
1044 384 : coord_map(i) = i + beglatxy - 2
1045 : else
1046 4224 : coord_map(i) = 0
1047 : end if
1048 : end do
1049 1536 : if (iam .lt. grid%npes_xy) then
1050 4608 : allocate(latvals(beglatxy:endlatxy))
1051 1536 : if (beglatxy == 1) then
1052 24 : latvals(1) = 0
1053 96 : latvals(2:endlatxy) = latdeg_st(1:endlatxy-1)
1054 : else
1055 1512 : i = beglatxy - 1 ! Stupid NAG 'error'
1056 7560 : latvals(beglatxy:endlatxy) = latdeg_st(i:endlatxy-1)
1057 : end if
1058 : else
1059 0 : allocate(latvals(0))
1060 : end if
1061 : slat_coord => horiz_coord_create('slat', '', (plat - 1), &
1062 : 'staggered latitude', 'degrees_north', beglatxy, endlatxy, latvals, &
1063 1536 : map=coord_map)
1064 1536 : deallocate(coord_map)
1065 1536 : deallocate(latvals)
1066 :
1067 : call cam_grid_register('fv_u_stagger', dyn_ustag_decomp, slat_coord, &
1068 1536 : lon_coord, grid_map, unstruct=.false.)
1069 : call cam_grid_attribute_register('fv_u_stagger', 'w_stag', &
1070 1536 : 'staggered latitude weights', 'slat', w_staggered)
1071 1536 : nullify(grid_map) ! Belongs to the grid
1072 :
1073 : ! Staggered grid for V_S
1074 : ! Make a map (need to do this because lat indices are distributed)
1075 4608 : allocate(grid_map(4, ((endlonxy - beglonxy + 1) * (endlatxy - beglatxy + 1))))
1076 1536 : ind = 0
1077 6144 : do i = beglatxy, endlatxy
1078 116736 : do j = beglonxy, endlonxy
1079 110592 : ind = ind + 1
1080 110592 : grid_map(1, ind) = j
1081 110592 : grid_map(2, ind) = i
1082 110592 : grid_map(3, ind) = j
1083 115200 : grid_map(4, ind) = i
1084 : end do
1085 : end do
1086 : ! Staggered longitude coordinate
1087 1536 : if (is_lon_distributed) then
1088 4608 : allocate(coord_map(endlonxy - beglonxy + 1))
1089 1536 : if (endlonxy >= beglonxy) then
1090 1536 : if (beglatxy == 1) then
1091 1776 : coord_map = (/ (i, i = beglonxy, endlonxy) /)
1092 : else
1093 37800 : coord_map = 0
1094 : end if
1095 : end if
1096 : slon_coord => horiz_coord_create('slon', '', plon, 'staggered longitude', &
1097 0 : 'degrees_east', beglonxy, endlonxy, londeg_st(beglonxy:endlonxy,1), &
1098 1536 : map=coord_map)
1099 1536 : deallocate(coord_map)
1100 : else
1101 : slon_coord => horiz_coord_create('slon', '', plon, 'staggered longitude', &
1102 0 : 'degrees_east', beglonxy, endlonxy, londeg_st(beglonxy:endlonxy,1))
1103 : end if
1104 : call cam_grid_register('fv_v_stagger', dyn_vstag_decomp, lat_coord, &
1105 1536 : slon_coord, grid_map, unstruct=.false.)
1106 1536 : nullify(grid_map) ! Belongs to the grid
1107 :
1108 : ! Zonal mean grid
1109 : ! Make a map
1110 4608 : allocate(grid_map(4, (endlatxy - beglatxy + 1)))
1111 1536 : ind = 0
1112 6144 : do i = beglatxy, endlatxy
1113 4608 : ind = ind + 1
1114 4608 : grid_map(1, ind) = 1
1115 4608 : grid_map(2, ind) = i
1116 6144 : if (beglonxy == 1) then
1117 384 : grid_map(3, ind) = 1
1118 384 : grid_map(4, ind) = i
1119 : else
1120 4224 : grid_map(3, ind) = 0
1121 4224 : grid_map(4, ind) = 0
1122 : end if
1123 : end do
1124 : ! We need a special, size-one "longigude" coordinate
1125 : ! NB: This is never a distributed coordinate so calc even on inactive PEs
1126 85231104 : zlon_bnds(1,1) = minval(londeg)
1127 85231104 : zlon_bnds(2,1) = maxval(londeg)
1128 1536 : allocate(latvals(1)) ! Really for a longitude
1129 1536 : latvals(1) = 0._r8
1130 : zlon_coord => horiz_coord_create('zlon', '', 1, 'longitude', &
1131 1536 : 'degrees_east', 1, 1, latvals(1:1), bnds=zlon_bnds)
1132 1536 : deallocate(latvals)
1133 : ! Zonal mean grid
1134 : call cam_grid_register('fv_centers_zonal', dyn_zonal_decomp, lat_coord, &
1135 1536 : zlon_coord, grid_map, unstruct=.false., zonal_grid=.true.)
1136 : ! Make sure 'gw' attribute shows up even if all variables are zonal mean
1137 1536 : call cam_grid_attribute_copy('fv_centers', 'fv_centers_zonal', 'gw')
1138 1536 : nullify(grid_map) ! Belongs to the grid
1139 :
1140 3072 : end subroutine define_cam_grids
1141 :
1142 : !========================================================================================
1143 :
1144 1536 : subroutine physgrid_copy_attributes_d(gridname, grid_attribute_names)
1145 1536 : use cam_grid_support, only: max_hcoordname_len
1146 :
1147 : ! Dummy arguments
1148 : character(len=max_hcoordname_len), intent(out) :: gridname
1149 : character(len=max_hcoordname_len), pointer, intent(out) :: grid_attribute_names(:)
1150 :
1151 1536 : gridname = 'fv_centers'
1152 1536 : allocate(grid_attribute_names(1))
1153 1536 : grid_attribute_names(1) = 'gw'
1154 :
1155 1536 : end subroutine physgrid_copy_attributes_d
1156 :
1157 : !========================================================================================
1158 :
1159 : end module dyn_grid
|
