LCOV - code coverage report
Current view: | top level - dynamics/se - restart_dynamics.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 276 | 433 | 63.7 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 7 | 11 | 63.6 %

              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       18432 :       do i = 1, ntrac
         162       16896 :          ierr = PIO_Def_Var(File, trim(cnst_name(i))//"_fvm", pio_double, &
         163       86016 :             (/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       18432 :       do m = 1, ntrac
         254      135696 :          do ie = 1, nelemd
         255    11184096 :             do k = 1, nlev
         256             :                ii = 1
         257    44312400 :                do j = 1, nc
         258   143629200 :                   do i = 1, nc
         259    99435600 :                      buf3d(ii,k,ie) = fvm(ie)%c(i,j,k,m)
         260   132580800 :                      ii = ii + 1
         261             :                   end do
         262             :                end do
         263             :             end do
         264             :          end do
         265       16896 :          call PIO_Setframe(file, c_fvm_desc(m), t_idx)
         266       18432 :          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        9216 :       do m = 1, ntrac
         668        8448 :          call infld(trim(cnst_name(m))//'_fvm', file, 'ncol_fvm', 'lev', &
         669             :                     1, nc**2, 1, nlev, 1, nelemd, var3d_fvm, readvar,    &
         670       16896 :                     gridname='FVM', timelevel=int(t_idx))
         671       68616 :          do ie = 1, nelemd
         672     5592048 :             do k = 1, nlev
         673             :                ii = 1
         674    22156200 :                do j = 1, nc
         675    71814600 :                   do i = 1, nc
         676    49717800 :                      fvm(ie)%c(i,j,k,m) = var3d_fvm(ii,k,ie)
         677    66290400 :                      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

Generated by: LCOV version 1.14