Line data Source code
1 : module restart_dynamics
2 :
3 : ! Write and read dynamics fields from the restart file. For exact restart
4 : ! it is necessary to write all element data, including duplicate columns,
5 : ! to the file. The namelist option, se_write_restart_unstruct, is
6 : ! available to write just the unique columns to the restart file using the
7 : ! same unstructured grid used by the history and initial files. This
8 : ! results in the introduction of a roundoff size difference on restart, but
9 : ! writes the fields in the unstructured grid format which is easier to
10 : ! modify if the user desires to introduce perturbations or other
11 : ! adjustments into the run. The restart file containing the unstructured
12 : ! grid format may also be used for an initial run.
13 :
14 : use shr_kind_mod, only: r8 => shr_kind_r8
15 : use spmd_utils, only: iam, masterproc
16 :
17 : use constituents, only: cnst_name
18 : use dyn_grid, only: timelevel, fvm, elem, edgebuf
19 : use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_init, write_restart_unstruct
20 : use hycoef, only: init_restart_hycoef, write_restart_hycoef, &
21 : hyai, hybi, ps0
22 : use ref_pres, only: ptop_ref
23 :
24 : use pio, only: pio_global, pio_unlimited, pio_offset_kind, pio_int, pio_double, &
25 : pio_seterrorhandling, pio_bcast_error, pio_noerr, &
26 : file_desc_t, var_desc_t, io_desc_t, &
27 : pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, &
28 : pio_def_dim, pio_def_var, &
29 : pio_enddef, &
30 : pio_initdecomp, pio_freedecomp, pio_setframe, &
31 : pio_put_att, pio_put_var, pio_write_darray, &
32 : pio_get_att, pio_get_var, pio_read_darray
33 :
34 : use cam_pio_utils, only: pio_subsystem, cam_pio_handle_error
35 : use cam_grid_support, only: cam_grid_header_info_t, cam_grid_id, cam_grid_write_attr, &
36 : cam_grid_write_var, cam_grid_get_decomp, cam_grid_dimensions, &
37 : max_hcoordname_len, cam_grid_get_dim_names
38 : use ncdio_atm, only: infld
39 :
40 : use infnan, only: isnan
41 : use cam_logfile, only: iulog
42 : use cam_abortutils, only: endrun
43 :
44 : use parallel_mod, only: par
45 : use thread_mod, only: horz_num_threads
46 : use dimensions_mod, only: np, npsq, ne, nlev, qsize, nelemd, nc, ntrac, use_cslam
47 : use dof_mod, only: UniquePoints
48 : use element_mod, only: element_t
49 : use se_dyn_time_mod, only: tstep, TimeLevel_Qdp
50 :
51 : use edge_mod, only: initEdgeBuffer, edgeVpack, edgeVunpack, FreeEdgeBuffer
52 : use edgetype_mod, only: EdgeBuffer_t
53 : use bndry_mod, only: bndry_exchange
54 :
55 : use fvm_control_volume_mod, only: fvm_struct
56 :
57 : implicit none
58 : private
59 : save
60 :
61 : public :: &
62 : init_restart_dynamics, &
63 : write_restart_dynamics, &
64 : read_restart_dynamics
65 :
66 : ! these variables are module data so they can be shared between the
67 : ! file definition and write phases
68 : type(var_desc_t) :: psdry_desc, udesc, vdesc, tdesc
69 : type(var_desc_t), allocatable :: qdesc_dp(:)
70 : type(var_desc_t) :: dp_fvm_desc
71 : type(var_desc_t), pointer :: c_fvm_desc(:)
72 :
73 : integer, private :: nelem_tot = -1 ! Correct total number of elements
74 :
75 : !=========================================================================================
76 : CONTAINS
77 : !=========================================================================================
78 :
79 2304 : subroutine init_nelem_tot()
80 : use spmd_utils, only: mpicom, MPI_INTEGER, MPI_SUM
81 :
82 : integer :: ierr
83 :
84 2304 : if (nelem_tot < 0) then
85 1536 : call MPI_Allreduce(nelemd, nelem_tot, 1, MPI_INTEGER, MPI_SUM, mpicom, ierr)
86 : end if
87 2304 : end subroutine init_nelem_tot
88 :
89 1536 : subroutine init_restart_dynamics(file, dyn_out)
90 :
91 : ! Define dimensions, variables, attributes for restart file.
92 :
93 : ! This is not really an "init" routine. It is called before
94 : ! write_restart_dynamics every time an restart is written.
95 :
96 : ! arguments
97 : type(file_desc_t), intent(inout) :: file
98 : type(dyn_export_t), intent(in) :: dyn_out
99 :
100 : ! local variables
101 : integer :: i
102 : integer :: vdimids(2)
103 : integer :: nlev_dimid
104 : integer :: ncol_dimid
105 : integer :: ncol_fvm_dimid
106 : integer :: time_dimid
107 :
108 : integer :: ierr, err_handling
109 :
110 : integer :: grid_id
111 1536 : type(cam_grid_header_info_t) :: info
112 :
113 : !----------------------------------------------------------------------------
114 :
115 1536 : call init_nelem_tot()
116 1536 : call init_restart_hycoef(file, vdimids)
117 1536 : nlev_dimid = vdimids(1)
118 :
119 1536 : call pio_seterrorhandling(File, pio_bcast_error, err_handling)
120 :
121 1536 : ierr = PIO_Def_Dim(File, 'time', PIO_UNLIMITED, time_dimid)
122 :
123 : ! GLL restart fields
124 :
125 : ! number of columns written to restart depends on whether all columns in the
126 : ! element structures are written, or just the unique columns (unstructured grid)
127 1536 : if (write_restart_unstruct) then
128 0 : grid_id = cam_grid_id('GLL')
129 0 : call cam_grid_write_attr(File, grid_id, info)
130 0 : ncol_dimid = info%get_hdimid(1)
131 : else
132 1536 : ierr = PIO_Def_Dim(File,'nenpnp', nelem_tot*np*np, ncol_dimid)
133 1536 : ierr = PIO_Put_Att(File, PIO_GLOBAL, 'ne', ne)
134 1536 : ierr = PIO_Put_Att(File, PIO_GLOBAL, 'np', np)
135 : end if
136 :
137 4608 : ierr = PIO_Def_Var(File, 'PSDRY', pio_double, (/ncol_dimid, time_dimid/), psdry_desc)
138 6144 : ierr = PIO_Def_Var(File, 'U', pio_double, (/ncol_dimid, nlev_dimid, time_dimid/), Udesc)
139 6144 : ierr = PIO_Def_Var(File, 'V', pio_double, (/ncol_dimid, nlev_dimid, time_dimid/), Vdesc)
140 6144 : ierr = PIO_Def_Var(File, 'T', pio_double, (/ncol_dimid, nlev_dimid, time_dimid/), Tdesc)
141 :
142 4608 : allocate(qdesc_dp(qsize))
143 :
144 10752 : do i=1,qsize
145 9216 : ierr = PIO_Def_Var(File,"dp"//trim(cnst_name(i)), pio_double, &
146 47616 : (/ncol_dimid, nlev_dimid, time_dimid/), Qdesc_dp(i))
147 : end do
148 :
149 : ! CSLAM restart fields
150 :
151 1536 : if (use_cslam) then
152 :
153 1536 : grid_id = cam_grid_id('FVM')
154 1536 : call cam_grid_write_attr(File, grid_id, info)
155 1536 : ncol_fvm_dimid = info%get_hdimid(1)
156 :
157 : ierr = PIO_Def_Var(File, 'dp_fvm', pio_double, &
158 6144 : (/ncol_fvm_dimid, nlev_dimid, time_dimid/), dp_fvm_desc)
159 :
160 4608 : allocate(c_fvm_desc(ntrac))
161 64512 : do i = 1, ntrac
162 62976 : ierr = PIO_Def_Var(File, trim(cnst_name(i))//"_fvm", pio_double, &
163 316416 : (/ncol_fvm_dimid, nlev_dimid, time_dimid/), c_fvm_desc(i))
164 : end do
165 :
166 : end if
167 :
168 1536 : call pio_seterrorhandling(File, err_handling)
169 :
170 1536 : end subroutine init_restart_dynamics
171 :
172 : !=========================================================================================
173 :
174 1536 : subroutine write_restart_dynamics(File, dyn_out)
175 : use control_mod, only: qsplit
176 : type(file_desc_t), intent(inout) :: File
177 : type(dyn_export_t), intent(in) :: dyn_out
178 :
179 : ! local variables
180 : integer(pio_offset_kind), parameter :: t_idx = 1
181 :
182 1536 : type(element_t), pointer :: elem(:)
183 1536 : type(fvm_struct), pointer :: fvm(:)
184 :
185 : integer :: tl, tlqdp
186 : integer :: i, ie, ii, j, k, m
187 : integer :: ierr
188 :
189 : integer :: grid_id
190 : integer :: grid_dimlens(2)
191 :
192 :
193 :
194 : integer :: array_lens(3)
195 : integer :: file_lens(2)
196 : type(io_desc_t), pointer :: iodesc3d_fvm
197 1536 : real(r8), allocatable :: buf3d(:,:,:)
198 :
199 :
200 :
201 : character(len=*), parameter :: sub = 'write_restart_dynamics'
202 : !----------------------------------------------------------------------------
203 :
204 1536 : call write_restart_hycoef(File)
205 :
206 1536 : tl = timelevel%n0
207 1536 : call TimeLevel_Qdp(timelevel, qsplit, tlQdp)
208 :
209 1536 : if (iam .lt. par%nprocs) then
210 1536 : elem => dyn_out%elem
211 1536 : fvm => dyn_out%fvm
212 : else
213 0 : allocate (elem(0), fvm(0))
214 : endif
215 :
216 : ! write fields on GLL grid
217 :
218 1536 : if (write_restart_unstruct) then
219 0 : call write_unstruct()
220 : else
221 1536 : call write_elem()
222 : end if
223 :
224 : ! write CSLAM fields
225 :
226 1536 : if (use_cslam) then
227 :
228 1536 : grid_id = cam_grid_id('FVM')
229 :
230 : ! write coords for FVM grid
231 1536 : call cam_grid_write_var(File, grid_id)
232 :
233 1536 : call cam_grid_dimensions(grid_id, grid_dimlens)
234 4608 : allocate(buf3d(nc*nc,nlev,nelemd))
235 6144 : array_lens = (/nc*nc, nlev, nelemd/)
236 4608 : file_lens = (/grid_dimlens(1), nlev/)
237 1536 : call cam_grid_get_decomp(grid_id, array_lens, file_lens, pio_double, iodesc3d_fvm)
238 :
239 12336 : do ie = 1, nelemd
240 1016736 : do k = 1, nlev
241 : ii = 1
242 4028400 : do j = 1, nc
243 13057200 : do i = 1, nc
244 9039600 : buf3d(ii,k,ie) = fvm(ie)%dp_fvm(i,j,k)
245 12052800 : ii = ii + 1
246 : end do
247 : end do
248 : end do
249 : end do
250 1536 : call PIO_Setframe(file, dp_fvm_desc, t_idx)
251 1536 : call PIO_Write_Darray(file, dp_fvm_desc, iodesc3d_fvm, buf3d, ierr)
252 :
253 64512 : do m = 1, ntrac
254 505776 : do ie = 1, nelemd
255 41686176 : do k = 1, nlev
256 : ii = 1
257 165164400 : do j = 1, nc
258 535345200 : do i = 1, nc
259 370623600 : buf3d(ii,k,ie) = fvm(ie)%c(i,j,k,m)
260 494164800 : ii = ii + 1
261 : end do
262 : end do
263 : end do
264 : end do
265 62976 : call PIO_Setframe(file, c_fvm_desc(m), t_idx)
266 64512 : call PIO_Write_Darray(file, c_fvm_desc(m), iodesc3d_fvm, buf3d, ierr)
267 : end do
268 :
269 1536 : deallocate(c_fvm_desc)
270 3072 : deallocate(buf3d)
271 : ! should this call be made on a pointer?
272 : !call pio_freedecomp(File, iodesc3d_fvm)
273 :
274 : end if
275 :
276 3072 : if (iam >= par%nprocs) then
277 0 : deallocate(elem, fvm)
278 : endif
279 :
280 : !-------------------------------------------------------------------------------
281 : contains
282 : !-------------------------------------------------------------------------------
283 :
284 1536 : subroutine write_elem()
285 :
286 : ! local variables
287 : integer :: i, ie, j, k
288 : integer :: ierr
289 1536 : integer, pointer :: ldof(:)
290 :
291 : type(io_desc_t) :: iodesc2d, iodesc3d
292 :
293 1536 : real(kind=r8), pointer :: var3d(:,:,:,:), var2d(:,:,:)
294 : !----------------------------------------------------------------------------
295 :
296 3072 : ldof => get_restart_decomp(elem, 1)
297 3072 : call PIO_InitDecomp(pio_subsystem, pio_double, (/nelem_tot*np*np/), ldof, iodesc2d)
298 1536 : deallocate(ldof)
299 :
300 1536 : ldof => get_restart_decomp(elem, nlev)
301 4608 : call PIO_InitDecomp(pio_subsystem, pio_double, (/nelem_tot*np*np,nlev/), ldof, iodesc3d)
302 1536 : deallocate(ldof)
303 :
304 4608 : allocate(var2d(np,np,nelemd))
305 6144 : allocate(var3d(np,np,nelemd,nlev))
306 :
307 : !$omp parallel do num_threads(horz_num_threads) private(ie, j, i)
308 12336 : do ie = 1, nelemd
309 55536 : do j = 1, np
310 226800 : do i = 1, np
311 216000 : var2d(i,j,ie) = elem(ie)%state%psdry(i,j)
312 : end do
313 : end do
314 : end do
315 1536 : call PIO_Setframe(File, psdry_desc, t_idx)
316 1536 : call PIO_Write_Darray(File, psdry_desc, iodesc2d, var2d, ierr)
317 :
318 : !$omp parallel do num_threads(horz_num_threads) private(ie, k, j, i)
319 12336 : do ie = 1, nelemd
320 1016736 : do k = 1, nlev
321 5032800 : do j = 1, np
322 21092400 : do i = 1, np
323 20088000 : var3d(i,j,ie,k) = elem(ie)%state%V(i,j,1,k,tl)
324 : end do
325 : end do
326 : end do
327 : end do
328 1536 : call PIO_Setframe(File, Udesc, t_idx)
329 1536 : call PIO_Write_Darray(File, Udesc, iodesc3d, var3d, ierr)
330 :
331 : !$omp parallel do num_threads(horz_num_threads) private(ie, k, j, i)
332 12336 : do ie = 1, nelemd
333 1016736 : do k = 1, nlev
334 5032800 : do j = 1, np
335 21092400 : do i = 1, np
336 20088000 : var3d(i,j,ie,k) = elem(ie)%state%V(i,j,2,k,tl)
337 : end do
338 : end do
339 : end do
340 : end do
341 1536 : call PIO_Setframe(File, Vdesc, t_idx)
342 1536 : call PIO_Write_Darray(File, Vdesc, iodesc3d, var3d, ierr)
343 :
344 : !$omp parallel do num_threads(horz_num_threads) private(ie, k, j, i)
345 12336 : do ie = 1, nelemd
346 1016736 : do k = 1, nlev
347 5032800 : do j = 1, np
348 21092400 : do i = 1, np
349 20088000 : var3d(i,j,ie,k) = elem(ie)%state%T(i,j,k,tl)
350 : end do
351 : end do
352 : end do
353 : end do
354 1536 : call PIO_Setframe(File, Tdesc, t_idx)
355 1536 : call PIO_Write_Darray(File, Tdesc, iodesc3d, var3d, ierr)
356 :
357 10752 : do m = 1, qsize
358 :
359 : !$omp parallel do num_threads(horz_num_threads) private(ie, k, j, i)
360 74016 : do ie = 1, nelemd
361 6100416 : do k = 1, nlev
362 30196800 : do j = 1, np
363 126554400 : do i = 1, np
364 120528000 : var3d(i,j,ie,k) = elem(ie)%state%Qdp(i,j,k,m,tlQdp)
365 : end do
366 : end do
367 : end do
368 : end do
369 9216 : call PIO_Setframe(File, Qdesc_dp(m), t_idx)
370 10752 : call PIO_Write_Darray(File, Qdesc_dp(m), iodesc3d, var3d, ierr)
371 :
372 : end do
373 :
374 1536 : deallocate(var2d)
375 1536 : deallocate(var3d)
376 1536 : deallocate(qdesc_dp)
377 :
378 1536 : call pio_freedecomp(File, iodesc2d)
379 1536 : call pio_freedecomp(File, iodesc3d)
380 :
381 7680 : end subroutine write_elem
382 :
383 : !-------------------------------------------------------------------------------
384 :
385 0 : subroutine write_unstruct()
386 :
387 : ! local variables
388 : integer :: i, ie, ii, j, k
389 : integer :: ierr
390 :
391 : integer :: array_lens_3d(3), array_lens_2d(2)
392 : integer :: file_lens_2d(2), file_lens_1d(1)
393 :
394 : type(io_desc_t), pointer :: iodesc
395 0 : real(r8), allocatable :: var2d(:,:), var3d(:,:,:)
396 : !----------------------------------------------------------------------------
397 :
398 0 : grid_id = cam_grid_id('GLL')
399 :
400 : ! write coordinate variables for unstructured GLL grid
401 0 : call cam_grid_write_var(File, grid_id)
402 :
403 : ! create map for distributed write
404 0 : call cam_grid_dimensions(grid_id, grid_dimlens)
405 :
406 : ! create map for distributed write of 2D fields
407 0 : array_lens_2d = (/npsq, nelemd/)
408 0 : file_lens_1d = (/grid_dimlens(1)/)
409 0 : call cam_grid_get_decomp(grid_id, array_lens_2d, file_lens_1d, pio_double, iodesc)
410 :
411 0 : allocate(var2d(npsq,nelemd))
412 :
413 0 : do ie = 1, nelemd
414 : ii = 1
415 0 : do j = 1, np
416 0 : do i = 1, np
417 0 : var2d(ii,ie) = elem(ie)%state%psdry(i,j)
418 0 : ii = ii + 1
419 : end do
420 : end do
421 : end do
422 0 : call PIO_Setframe(File, psdry_desc, t_idx)
423 0 : call PIO_Write_Darray(File, psdry_desc, iodesc, var2d, ierr)
424 :
425 0 : nullify(iodesc)
426 0 : deallocate(var2d)
427 :
428 : ! create map for distributed write of 3D fields
429 0 : array_lens_3d = (/npsq, nlev, nelemd/)
430 0 : file_lens_2d = (/grid_dimlens(1), nlev/)
431 0 : call cam_grid_get_decomp(grid_id, array_lens_3d, file_lens_2d, pio_double, iodesc)
432 :
433 0 : allocate(var3d(npsq,nlev,nelemd))
434 :
435 0 : do ie = 1, nelemd
436 0 : do k = 1, nlev
437 : ii = 1
438 0 : do j = 1, np
439 0 : do i = 1, np
440 0 : var3d(ii,k,ie) = elem(ie)%state%V(i,j,1,k,tl)
441 0 : ii = ii + 1
442 : end do
443 : end do
444 : end do
445 : end do
446 0 : call PIO_Setframe(File, Udesc, t_idx)
447 0 : call PIO_Write_Darray(File, Udesc, iodesc, var3d, ierr)
448 :
449 0 : do ie = 1, nelemd
450 0 : do k = 1, nlev
451 : ii = 1
452 0 : do j = 1, np
453 0 : do i = 1, np
454 0 : var3d(ii,k,ie) = elem(ie)%state%V(i,j,2,k,tl)
455 0 : ii = ii + 1
456 : end do
457 : end do
458 : end do
459 : end do
460 0 : call PIO_Setframe(File, Vdesc, t_idx)
461 0 : call PIO_Write_Darray(File, Vdesc, iodesc, var3d, ierr)
462 :
463 0 : do ie = 1, nelemd
464 0 : do k = 1, nlev
465 : ii = 1
466 0 : do j = 1, np
467 0 : do i = 1, np
468 0 : var3d(ii,k,ie) = elem(ie)%state%T(i,j,k,tl)
469 0 : ii = ii + 1
470 : end do
471 : end do
472 : end do
473 : end do
474 0 : call PIO_Setframe(File, Tdesc, t_idx)
475 0 : call PIO_Write_Darray(File, Tdesc, iodesc, var3d, ierr)
476 :
477 0 : do m = 1, qsize
478 :
479 : !$omp parallel do num_threads(horz_num_threads) private(ie, k, j, i)
480 0 : do ie = 1, nelemd
481 0 : do k = 1, nlev
482 : ii = 1
483 0 : do j = 1, np
484 0 : do i = 1, np
485 0 : var3d(ii,k,ie) = elem(ie)%state%Qdp(i,j,k,m,tlQdp)
486 0 : ii = ii + 1
487 : end do
488 : end do
489 : end do
490 : end do
491 0 : call PIO_Setframe(File, Qdesc_dp(m), t_idx)
492 0 : call PIO_Write_Darray(File, Qdesc_dp(m), iodesc, var3d, ierr)
493 :
494 : end do
495 :
496 0 : deallocate(var3d)
497 0 : deallocate(qdesc_dp)
498 :
499 0 : end subroutine write_unstruct
500 :
501 : !-------------------------------------------------------------------------------
502 :
503 : end subroutine write_restart_dynamics
504 :
505 : !=========================================================================================
506 :
507 768 : subroutine read_restart_dynamics(File, dyn_in, dyn_out)
508 : use control_mod, only: qsplit
509 : ! arguments
510 : type(File_desc_t), intent(inout) :: File
511 : type(dyn_import_t), intent(out) :: dyn_in
512 : type(dyn_export_t), intent(out) :: dyn_out
513 :
514 : ! local variables
515 : integer(pio_offset_kind), parameter :: t_idx = 1
516 :
517 : integer :: tl, tlQdp
518 : integer :: i, ie, ii, k, m, j
519 : integer :: ierr, err_handling
520 : integer :: fne, fnp, fnlev, fnc
521 : integer :: hdim_len, ncols_fvm
522 :
523 : integer :: nlev_dimid
524 : integer :: ncol_dimid
525 : integer :: ncol_fvm_dimid
526 :
527 : type(var_desc_t) :: udesc
528 : type(var_desc_t) :: vdesc
529 : type(var_desc_t) :: tdesc
530 : type(var_desc_t) :: psdry_desc
531 768 : type(var_desc_t), allocatable :: qdesc_dp(:)
532 :
533 : integer :: grid_id
534 : integer :: grid_dimlens(2)
535 : character(len=max_hcoordname_len) :: dimname1, dimname2
536 :
537 768 : real(r8), allocatable :: var3d_fvm(:,:,:)
538 :
539 : logical :: readvar
540 :
541 : character(len=*), parameter :: sub = 'read_restart_dynamics'
542 : !----------------------------------------------------------------------------
543 :
544 : ! Note1: the hybrid coefficients are read from the same location as for an
545 : ! initial run (e.g., dyn_grid_init).
546 :
547 : ! Note2: the dyn_in and dyn_out objects are not associated with the elem and fvm
548 : ! objects until dyn_init is called. Until the restart is better integrated
549 : ! into dyn_init we just access elem and fvm directly from the dyn_grid
550 : ! module.
551 :
552 768 : tl = timelevel%n0
553 768 : call TimeLevel_Qdp(timelevel, qsplit, tlQdp)
554 768 : call init_nelem_tot()
555 :
556 768 : call pio_seterrorhandling(File, pio_bcast_error, err_handling)
557 :
558 : ! some checks that the restart contains the same grid as the running model.
559 :
560 768 : ierr = PIO_Get_Att(File, PIO_GLOBAL, 'ne', fne)
561 768 : ierr = PIO_Get_Att(File, PIO_GLOBAL, 'np', fnp)
562 768 : if (ne /= fne .or. np /= fnp) then
563 0 : write(iulog,*) 'Restart file np or ne does not match model. np (file, model):', &
564 0 : fnp, np, ' ne (file, model) ', fne, ne
565 0 : call endrun(sub//': Restart file np or ne does not match model.')
566 : end if
567 :
568 768 : ierr = PIO_Inq_DimID(File, 'lev', nlev_dimid)
569 768 : ierr = PIO_Inq_dimlen(File, nlev_dimid, fnlev)
570 768 : if (nlev /= fnlev) then
571 0 : write(iulog,*) 'Restart file nlev does not match model. nlev (file, namelist):', &
572 0 : fnlev, nlev
573 0 : call endrun(sub//': Restart file nlev does not match model.')
574 : end if
575 :
576 : ! variable descriptors of required dynamics fields
577 768 : ierr = PIO_Inq_varid(File, 'U', udesc)
578 768 : call cam_pio_handle_error(ierr, sub//': cannot find U')
579 768 : ierr = PIO_Inq_varid(File, 'V', Vdesc)
580 768 : call cam_pio_handle_error(ierr, sub//': cannot find V')
581 768 : ierr = PIO_Inq_varid(File, 'T', tdesc)
582 768 : call cam_pio_handle_error(ierr, sub//': cannot find T')
583 768 : ierr = PIO_Inq_varid(File, 'PSDRY', psdry_desc)
584 768 : call cam_pio_handle_error(ierr, sub//': cannot find PSDRY')
585 2304 : allocate(qdesc_dp(qsize))
586 5376 : do m = 1, qsize
587 4608 : ierr = PIO_Inq_varid(File, "dp"//trim(cnst_name(m)), Qdesc_dp(m))
588 5376 : call cam_pio_handle_error(ierr, sub//': cannot find dp'//trim(cnst_name(m)))
589 : end do
590 :
591 : ! check whether the restart fields on the GLL grid contain unique columns
592 : ! or the element structure (nenpnp = nelem_tot*np*np columns)
593 :
594 768 : ierr = PIO_Inq_DimID(File, 'nenpnp', ncol_dimid)
595 768 : if (ierr == pio_noerr) then
596 :
597 768 : call read_elem()
598 :
599 : else
600 :
601 0 : call read_unstruct()
602 :
603 : end if
604 :
605 768 : deallocate(qdesc_dp)
606 :
607 : ! recompute dp3d from psdry
608 6168 : do ie = 1, nelemd
609 508368 : do k = 1, nlev
610 1506600 : elem(ie)%state%dp3d(:,:,k,tl) = ((hyai(k+1) - hyai(k))*ps0) + &
611 12058200 : ((hybi(k+1) - hybi(k))*elem(ie)%state%psdry(:,:))
612 : end do
613 : end do
614 :
615 : ! Seems like this initialization should be done somewhere else.
616 6168 : do ie = 1, nelemd
617 21600000 : elem(ie)%derived%fM = 0._r8
618 10551600 : elem(ie)%derived%fT = 0._r8
619 63315768 : elem(ie)%derived%fQ = 0._r8
620 : end do
621 :
622 : ! read cslam fields
623 :
624 768 : if (use_cslam) then
625 :
626 : ! Checks that file and model dimensions agree.
627 :
628 768 : ierr = PIO_Get_Att(File, PIO_GLOBAL, 'nc', fnc)
629 768 : if (nc /= fnc) then
630 0 : write(iulog,*) 'Restart file nc does not match model. nc (file, model):',fnc,nc,&
631 0 : ' ne (file, model) ', fne, ne
632 0 : call endrun(sub//': Restart file nc does not match model.')
633 : end if
634 :
635 768 : ierr = PIO_Inq_DimID(File, 'ncol_fvm', ncol_fvm_dimid)
636 768 : call cam_pio_handle_error(ierr, sub//': cannot find ncol_fvm')
637 768 : ierr = PIO_Inq_dimlen(File, ncol_fvm_dimid, ncols_fvm)
638 :
639 768 : grid_id = cam_grid_id('FVM')
640 768 : call cam_grid_dimensions(grid_id, grid_dimlens)
641 :
642 768 : if (ncols_fvm /= grid_dimlens(1)) then
643 0 : write(iulog,*) 'Restart file ncol_fvm does not match model. ncols_fvm (file, model):',&
644 0 : ncols_fvm, grid_dimlens(1)
645 0 : call endrun(sub//': Restart file ncols_fvm does not match model.')
646 : end if
647 :
648 2304 : allocate(var3d_fvm(nc*nc,nlev,nelemd))
649 5028168 : var3d_fvm = 0._r8
650 :
651 : ! dp_fvm
652 : call infld('dp_fvm', file, 'ncol_fvm', 'lev', 1, nc**2, 1, nlev, 1, nelemd, &
653 768 : var3d_fvm, readvar, gridname='FVM', timelevel=int(t_idx))
654 6168 : do ie = 1, nelemd
655 508368 : do k = 1, nlev
656 : ii = 1
657 2014200 : do j = 1, nc
658 6528600 : do i = 1, nc
659 4519800 : fvm(ie)%dp_fvm(i,j,k) = var3d_fvm(ii,k,ie)
660 6026400 : ii = ii + 1
661 : end do
662 : end do
663 : end do
664 : end do
665 :
666 : ! tracers
667 32256 : do m = 1, ntrac
668 31488 : call infld(trim(cnst_name(m))//'_fvm', file, 'ncol_fvm', 'lev', &
669 : 1, nc**2, 1, nlev, 1, nelemd, var3d_fvm, readvar, &
670 62976 : gridname='FVM', timelevel=int(t_idx))
671 253656 : do ie = 1, nelemd
672 20843088 : do k = 1, nlev
673 : ii = 1
674 82582200 : do j = 1, nc
675 267672600 : do i = 1, nc
676 185311800 : fvm(ie)%c(i,j,k,m) = var3d_fvm(ii,k,ie)
677 247082400 : ii = ii + 1
678 : end do
679 : end do
680 : end do
681 : end do
682 : end do
683 :
684 : ! compute dry surface pressure (a derived quantity)
685 6168 : do ie = 1, nelemd
686 22368 : do j = 1, nc
687 70200 : do i = 1, nc
688 4584600 : fvm(ie)%psc(i,j) = sum(fvm(ie)%dp_fvm(i,j,:)) + ptop_ref
689 : end do
690 : end do
691 : end do
692 :
693 : end if
694 :
695 768 : call pio_seterrorhandling(File, err_handling)
696 :
697 768 : call dyn_init(dyn_in, dyn_out)
698 :
699 : !-------------------------------------------------------------------------------
700 : contains
701 : !-------------------------------------------------------------------------------
702 :
703 768 : subroutine read_elem()
704 :
705 : ! local variables
706 : integer :: ierr
707 : integer :: ncol
708 : integer :: i, ie, ii, j, k, m
709 :
710 768 : integer, pointer :: ldof(:)
711 :
712 : type(io_desc_t) :: iodesc2d, iodesc3d
713 768 : real(r8), allocatable :: var3d(:), var2d(:)
714 :
715 : character(len=*), parameter :: sub='read_elem'
716 : !----------------------------------------------------------------------------
717 :
718 1536 : ierr = PIO_Inq_dimlen(File, ncol_dimid, ncol)
719 768 : call cam_pio_handle_error(ierr, sub//': reading nenpnp')
720 :
721 768 : ldof => get_restart_decomp(elem, 1)
722 1536 : call PIO_InitDecomp(pio_subsystem, pio_double, (/ncol/), ldof, iodesc2d)
723 768 : deallocate(ldof)
724 :
725 768 : ldof => get_restart_decomp(elem, nlev)
726 2304 : call PIO_InitDecomp(pio_subsystem, pio_double, (/ncol,nlev/), ldof, iodesc3d)
727 768 : deallocate(ldof)
728 :
729 2304 : allocate(var2d(nelemd*np*np), stat=ierr)
730 768 : if (ierr/=0) call endrun( sub//': not able to allocate var2d' )
731 87168 : var2d = 0._r8
732 :
733 768 : call pio_setframe(File, psdry_desc, t_idx)
734 768 : call pio_read_darray(File, psdry_desc, iodesc2d, var2d, ierr)
735 768 : call cam_pio_handle_error(ierr, sub//': reading PSDRY')
736 768 : ii = 0
737 6168 : do ie = 1, nelemd
738 27768 : do j = 1, np
739 113400 : do i = 1, np
740 86400 : ii = ii + 1
741 108000 : elem(ie)%state%psdry(i,j) = var2d(ii)
742 : end do
743 : end do
744 : end do
745 :
746 768 : deallocate(var2d)
747 2304 : allocate(var3d(nelemd*np*np*nlev), stat=ierr)
748 768 : if (ierr/=0) call endrun( sub//': not able to allocate var3d' )
749 8035968 : var3d = 0._r8
750 :
751 768 : call pio_setframe(File, udesc, t_idx)
752 768 : call pio_read_darray(File, udesc, iodesc3d, var3d, ierr)
753 768 : call cam_pio_handle_error(ierr, sub//': reading U')
754 768 : ii = 0
755 72192 : do k = 1, nlev
756 574392 : do ie = 1, nelemd
757 2582424 : do j = 1, np
758 10546200 : do i = 1, np
759 8035200 : ii = ii + 1
760 10044000 : elem(ie)%state%v(i,j,1,k,tl) = var3d(ii)
761 : end do
762 : end do
763 : end do
764 : end do
765 :
766 768 : call pio_setframe(File, vdesc, t_idx)
767 768 : call pio_read_darray(File, vdesc, iodesc3d, var3d, ierr)
768 768 : call cam_pio_handle_error(ierr, sub//': reading V')
769 768 : ii = 0
770 72192 : do k = 1, nlev
771 574392 : do ie = 1, nelemd
772 2582424 : do j = 1, np
773 10546200 : do i = 1, np
774 8035200 : ii = ii + 1
775 10044000 : elem(ie)%state%v(i,j,2,k,tl) = var3d(ii)
776 : end do
777 : end do
778 : end do
779 : end do
780 :
781 768 : call pio_setframe(File, tdesc, t_idx)
782 768 : call pio_read_darray(File, tdesc, iodesc3d, var3d, ierr)
783 768 : call cam_pio_handle_error(ierr, sub//': reading T')
784 768 : ii = 0
785 72192 : do k = 1, nlev
786 574392 : do ie = 1, nelemd
787 2582424 : do j = 1, np
788 10546200 : do i = 1, np
789 8035200 : ii = ii + 1
790 10044000 : elem(ie)%state%T(i,j,k,tl) = var3d(ii)
791 : end do
792 : end do
793 : end do
794 : end do
795 :
796 5376 : do m = 1, qsize
797 4608 : call pio_setframe(File, qdesc_dp(m), t_idx)
798 4608 : call pio_read_darray(File, qdesc_dp(m), iodesc3d, var3d, ierr)
799 4608 : call cam_pio_handle_error(ierr, sub//': reading dp'//trim(cnst_name(m)))
800 4608 : ii = 0
801 438528 : do k = 1, nlev
802 3446352 : do ie = 1, nelemd
803 15494544 : do j = 1, np
804 63277200 : do i = 1, np
805 48211200 : ii = ii + 1
806 60264000 : elem(ie)%state%Qdp(i,j,k,m,tlQdp) = var3d(ii)
807 : end do
808 : end do
809 : end do
810 : end do
811 : end do
812 :
813 768 : deallocate(var3d)
814 :
815 768 : if (masterproc) write(iulog,*) sub//': completed successfully'
816 :
817 4608 : end subroutine read_elem
818 :
819 : !-------------------------------------------------------------------------------
820 :
821 0 : subroutine read_unstruct()
822 :
823 : ! local variables
824 : integer :: grid_id
825 :
826 : integer :: i, ie, ii, j, kptr, m
827 :
828 0 : real(r8), allocatable :: dbuf2(:,:) ! (npsq,nelemd)
829 0 : real(r8), allocatable :: dbuf3(:,:,:) ! (npsq,nlev,nelemd)
830 :
831 0 : type(EdgeBuffer_t) :: edge
832 :
833 : character(len=*), parameter :: sub='read_unstruct'
834 : !----------------------------------------------------------------------------
835 :
836 : ! The name of the unstructured grid dimension is in the grid object
837 : ! since the coordinate date was written to the restart file using
838 : ! that object.
839 0 : grid_id = cam_grid_id('GLL')
840 0 : call cam_grid_get_dim_names(grid_id, dimname1, dimname2)
841 :
842 0 : allocate(dbuf2(npsq,nelemd))
843 :
844 0 : call read_2d('PSDRY', dbuf2)
845 0 : do ie = 1, nelemd
846 : ii = 1
847 0 : do j = 1, np
848 0 : do i = 1, np
849 0 : elem(ie)%state%psdry(i,j) = dbuf2(ii,ie)
850 0 : ii = ii + 1
851 : end do
852 : end do
853 : end do
854 :
855 0 : deallocate(dbuf2)
856 :
857 0 : allocate(dbuf3(npsq,nlev,nelemd))
858 :
859 0 : call read_3d('U', dbuf3)
860 0 : do ie = 1, nelemd
861 : ii = 1
862 0 : do j = 1, np
863 0 : do i = 1, np
864 0 : elem(ie)%state%v(i,j,1,:,tl) = dbuf3(ii,:,ie)
865 0 : ii = ii + 1
866 : end do
867 : end do
868 : end do
869 :
870 0 : call read_3d('V', dbuf3)
871 0 : do ie = 1, nelemd
872 : ii = 1
873 0 : do j = 1, np
874 0 : do i = 1, np
875 0 : elem(ie)%state%v(i,j,2,:,tl) = dbuf3(ii,:,ie)
876 0 : ii = ii + 1
877 : end do
878 : end do
879 : end do
880 :
881 0 : call read_3d('T', dbuf3)
882 0 : do ie = 1, nelemd
883 : ii = 1
884 0 : do j = 1, np
885 0 : do i = 1, np
886 0 : elem(ie)%state%T(i,j,:,tl) = dbuf3(ii,:,ie)
887 0 : ii = ii + 1
888 : end do
889 : end do
890 : end do
891 :
892 0 : do m = 1, qsize
893 :
894 0 : call read_3d('dp'//trim(cnst_name(m)), dbuf3)
895 0 : do ie = 1, nelemd
896 : ii = 1
897 0 : do j = 1, np
898 0 : do i = 1, np
899 0 : elem(ie)%state%Qdp(i,j,:,m,tlQdp) = dbuf3(ii,:,ie)
900 0 : ii = ii + 1
901 : end do
902 : end do
903 : end do
904 :
905 : end do
906 :
907 0 : deallocate(dbuf3)
908 :
909 : ! boundary exchange
910 0 : if (iam < par%nprocs) then
911 0 : call initEdgeBuffer(par, edge, elem, (3+qsize)*nlev + 1 )
912 : end if
913 0 : do ie = 1, nelemd
914 0 : kptr = 0
915 0 : call edgeVpack(edge, elem(ie)%state%psdry(:,:), 1, kptr, ie)
916 0 : kptr = kptr + 1
917 0 : call edgeVpack(edge, elem(ie)%state%v(:,:,:,:,tl), 2*nlev, kptr, ie)
918 0 : kptr = kptr + (2 * nlev)
919 0 : call edgeVpack(edge, elem(ie)%state%T(:,:,:,tl), nlev, kptr, ie)
920 0 : kptr = kptr + nlev
921 0 : call edgeVpack(edge, elem(ie)%state%Qdp(:,:,:,:,tlQdp), nlev*qsize, kptr, ie)
922 : end do
923 0 : if (iam < par%nprocs) then
924 0 : call bndry_exchange(par, edge, location='read_restart_dynamics::read_ustruct')
925 : end if
926 0 : do ie = 1, nelemd
927 0 : kptr = 0
928 0 : call edgeVunpack(edge, elem(ie)%state%psdry(:,:), 1, kptr, ie)
929 0 : kptr = kptr + 1
930 0 : call edgeVunpack(edge, elem(ie)%state%v(:,:,:,:,tl), 2*nlev, kptr, ie)
931 0 : kptr = kptr + (2 * nlev)
932 0 : call edgeVunpack(edge, elem(ie)%state%T(:,:,:,tl), nlev, kptr, ie)
933 0 : kptr = kptr + nlev
934 0 : call edgeVunpack(edge, elem(ie)%state%Qdp(:,:,:,:,tlQdp), nlev*qsize, kptr, ie)
935 : end do
936 :
937 0 : if (iam < par%nprocs) then
938 0 : call FreeEdgeBuffer(edge)
939 : end if
940 :
941 0 : end subroutine read_unstruct
942 :
943 : !-------------------------------------------------------------------------------
944 :
945 0 : subroutine read_2d(fieldname, buffer)
946 :
947 : character(len=*), intent(in) :: fieldname
948 : real(r8), intent(inout) :: buffer(:, :)
949 :
950 : logical :: found
951 : !----------------------------------------------------------------------------
952 :
953 0 : buffer = 0.0_r8
954 : call infld(trim(fieldname), file, dimname1, 1, npsq, 1, nelemd, buffer, &
955 0 : found, gridname='GLL', timelevel=int(t_idx))
956 0 : if (.not. found) then
957 : call endrun('read_restart_dynamics: read_unstruct: read_2d: Could not find ' // &
958 0 : trim(fieldname))
959 : end if
960 :
961 : ! This code allows use of compiler option to set uninitialized values
962 : ! to NaN. In that case infld can return NaNs where the element GLL points
963 : ! are not "unique columns"
964 0 : where (isnan(buffer)) buffer = 0.0_r8
965 :
966 0 : end subroutine read_2d
967 :
968 : !-------------------------------------------------------------------------------
969 :
970 0 : subroutine read_3d(fieldname, buffer)
971 :
972 : character(len=*), intent(in) :: fieldname
973 : real(r8), intent(inout) :: buffer(:,:,:)
974 :
975 : logical :: found
976 : !----------------------------------------------------------------------------
977 :
978 0 : buffer = 0.0_r8
979 : call infld(trim(fieldname), file, dimname1, 'lev', 1, npsq, 1, nlev, &
980 0 : 1, nelemd, buffer, found, gridname='GLL', timelevel=int(t_idx))
981 0 : if (.not. found) then
982 : call endrun('read_restart_dynamics: read_unstruct: read_3d: Could not find ' // &
983 0 : trim(fieldname))
984 : end if
985 :
986 : ! This code allows use of compiler option to set uninitialized values
987 : ! to NaN. In that case infld can return NaNs where the element GLL points
988 : ! are not "unique columns"
989 0 : where (isnan(buffer)) buffer = 0.0_r8
990 :
991 0 : end subroutine read_3d
992 :
993 : !-------------------------------------------------------------------------------
994 : end subroutine read_restart_dynamics
995 :
996 : !=========================================================================================
997 : ! Private
998 : !=========================================================================================
999 :
1000 4608 : function get_restart_decomp(elem, lev) result(ldof)
1001 :
1002 : ! Get the integer mapping of a variable in the dynamics decomp in memory.
1003 : ! The canonical ordering is as on the file. A 0 value indicates that the
1004 : ! variable is not on the file (eg halo or boundary values)
1005 :
1006 : type(element_t), intent(in) :: elem(:)
1007 : integer, intent(in) :: lev
1008 : integer, pointer :: ldof(:)
1009 :
1010 : integer :: i, j, k, ie
1011 : !----------------------------------------------------------------------------
1012 :
1013 13824 : allocate(ldof(nelemd*np*np*lev))
1014 :
1015 4608 : j = 1
1016 221184 : do k = 1, lev
1017 1743984 : do ie = 1, nelemd
1018 26104176 : do i = 1, np*np
1019 24364800 : ldof(j) = (elem(ie)%GlobalID-1)*np*np + (k-1)*nelem_tot*np*np + i
1020 25887600 : j = j + 1
1021 : end do
1022 : end do
1023 : end do
1024 :
1025 4608 : end function get_restart_decomp
1026 :
1027 : !=========================================================================================
1028 :
1029 : function get_restart_decomp_fvm(elem, lev) result(ldof)
1030 :
1031 : type(element_t), intent(in) :: elem(:)
1032 : integer, intent(in) :: lev
1033 : integer, pointer :: ldof(:)
1034 :
1035 : integer :: i, j, k, ie
1036 : !----------------------------------------------------------------------------
1037 :
1038 : allocate(ldof(nelemd*nc*nc*lev))
1039 :
1040 : j = 1
1041 : do k = 1, lev
1042 : do ie = 1, nelemd
1043 : do i = 1, nc*nc
1044 : ldof(j) = (elem(ie)%GlobalID-1)*nc*nc + (k-1)*nelem_tot*nc*nc + i
1045 : j = j + 1
1046 : end do
1047 : end do
1048 : end do
1049 :
1050 : end function get_restart_decomp_fvm
1051 :
1052 : !=========================================================================================
1053 :
1054 : end module restart_dynamics
|
