LCOV - code coverage report
Current view: | top level - control - ncdio_atm.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 94 | 291 | 32.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 4 | 4 | 100.0 %

              Line data    Source code
           1             : module ncdio_atm
           2             :
           3             :   !-----------------------------------------------------------------------
           4             :   !BOP
           5             :   !
           6             :   ! !MODULE: ncdio_atm
           7             :   !
           8             :   ! !DESCRIPTION:
           9             :   ! Generic interfaces to write fields to PIO files
          10             :   !
          11             :   ! !USES:
          12             :
          13             :   use pio, only: pio_offset_kind, file_desc_t, var_desc_t, pio_double,        &
          14             :        pio_inq_dimid, pio_max_var_dims, io_desc_t, pio_setframe
          15             :   use shr_kind_mod,   only: r8 => shr_kind_r8
          16             :   use shr_sys_mod,    only: shr_sys_flush      ! Standardized system subroutines
          17             :   use shr_scam_mod,   only: shr_scam_getCloseLatLon  ! Standardized system subroutines
          18             :   use spmd_utils,     only: masterproc
          19             :   use cam_abortutils, only: endrun
          20             :   use scamMod,        only: scmlat,scmlon,single_column
          21             :   use cam_logfile,    only: iulog
          22             :   use string_utils,   only: to_lower
          23             :   use cam_grid_support, only: cam_grid_check, cam_grid_get_decomp, cam_grid_id, &
          24             :                               cam_grid_dimensions, cam_grid_get_latvals, cam_grid_get_lonvals, &
          25             :                               max_hcoordname_len
          26             :   !
          27             :   ! !PUBLIC TYPES:
          28             :   implicit none
          29             :
          30             :   PRIVATE
          31             :
          32             :   save
          33             :
          34             :   logical :: debug = .false.
          35             :
          36             :   !
          37             :   !EOP
          38             :   !
          39             :   interface infld
          40             :      module procedure infld_real_1d_2d
          41             :      module procedure infld_real_2d_2d
          42             :      module procedure infld_real_2d_3d
          43             :      module procedure infld_real_3d_3d
          44             :   end interface
          45             :
          46             :   public :: infld
          47             :
          48             :   !-----------------------------------------------------------------------
          49             :
          50             : contains
          51             :
          52             :   !-----------------------------------------------------------------------
          53             :   !BOP
          54             :   !
          55             :   ! !IROUTINE: infld_real_1d_2d
          56             :   !
          57             :   ! !INTERFACE:
          58      204288 :   subroutine infld_real_1d_2d(varname, ncid, dimname1,                        &
          59      146688 :        dim1b, dim1e, dim2b, dim2e, field, readvar, gridname, timelevel,       &
          60             :        fillvalue)
          61             :     !
          62             :     ! !DESCRIPTION:
          63             :     ! Netcdf I/O of initial real field from netCDF file
          64             :     ! Read a 1-D field (or slice) into a 2-D variable
          65             :     !
          66             :     ! !USES
          67             :     !
          68             :
          69             :     use pio,              only: pio_read_darray, pio_setdebuglevel
          70             :     use pio,              only: PIO_MAX_NAME, pio_inq_dimname
          71             :     use cam_pio_utils,    only: cam_pio_check_var, cam_pio_inq_var_fill
          72             :
          73             :     !
          74             :     ! !ARGUMENTS:
          75             :     implicit none
          76             :     character(len=*),  intent(in)     :: varname  ! variable name
          77             :     type(file_desc_t), intent(inout)  :: ncid     ! input unit
          78             :     character(len=*),  intent(in)     :: dimname1 ! name of 1st array dimensions of field on file (array order)
          79             :     integer,           intent(in)     :: dim1b    ! start of first  dimension of array to be returned
          80             :     integer,           intent(in)     :: dim1e    ! end   of first  dimension of array to be returned
          81             :     integer,           intent(in)     :: dim2b    ! start of second dimension of array to be returned
          82             :     integer,           intent(in)     :: dim2e    ! end   of second dimension of array to be returned
          83             :     real(r8), target,  intent(out)    :: field(dim1b:dim1e,dim2b:dim2e) ! array to be returned (decomposed or global)
          84             :     logical,           intent(out)    :: readvar  ! true => variable is on initial dataset
          85             :     character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
          86             :     integer, optional, intent(in)     :: timelevel
          87             :     real(r8), optional, intent(out)   :: fillvalue
          88             :     !
          89             :     !EOP
          90             :     !
          91             :     ! !LOCAL VARIABLES:
          92             :     type(io_desc_t), pointer  :: iodesc
          93             :     integer                   :: grid_id   ! grid ID for data mapping
          94             :     integer                   :: j         ! index
          95             :     integer                   :: ierr      ! error status
          96             :     type(var_desc_t)          :: varid     ! variable id
          97             :     integer                   :: no_fill
          98             :     integer                   :: arraydimsize(2) ! field dimension lengths
          99             :
         100             :     integer                   :: ndims ! number of dimensions
         101             :     integer                   :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
         102             :     integer                   :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
         103             :     integer                   :: grid_dimlens(2)
         104             :
         105             :     character(len=PIO_MAX_NAME) :: tmpname
         106             :     character(len=128)        :: errormsg
         107             :
         108             :     logical                   :: readvar_tmp ! if true, variable is on tape
         109             :     character(len=*), parameter :: subname='INFLD_REAL_1D_2D' ! subroutine name
         110             :     character(len=max_hcoordname_len) :: vargridname ! Name of variable's grid
         111             :
         112             :     !
         113             :     !-----------------------------------------------------------------------
         114             :     !
         115             :     !    call pio_setdebuglevel(3)
         116             :
         117      146688 :     nullify(iodesc)
         118             :
         119             :     !
         120             :     ! Error conditions
         121             :     !
         122      146688 :     if (present(gridname)) then
         123      146688 :        vargridname=trim(gridname)
         124             :     else
         125           0 :        vargridname='physgrid'
         126             :     end if
         127             :
         128      146688 :     if (single_column .and. vargridname=='physgrid') then
         129           0 :        vargridname='physgrid_scm'
         130             :     end if
         131             :
         132      146688 :     grid_id = cam_grid_id(trim(vargridname))
         133             :
         134      146688 :     if (.not. cam_grid_check(grid_id)) then
         135           0 :       if(masterproc) then
         136           0 :          write(errormsg, *)': invalid gridname, "',trim(vargridname),'", specified for field ',trim(varname)
         137             :       end if
         138           0 :       call endrun(trim(subname)//errormsg)
         139             :     end if
         140             :
         141      146688 :     if (debug .and. masterproc) then
         142           0 :        write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(vargridname)
         143           0 :        call shr_sys_flush(iulog)
         144             :     end if
         145             :
         146             :     ! Get the number of columns in the global grid.
         147      146688 :     call cam_grid_dimensions(grid_id, grid_dimlens)
         148             :     !
         149             :     ! Read netCDF file
         150             :     !
         151             :     !
         152             :     ! Check if field is on file; get netCDF variable id
         153             :     !
         154      146688 :     call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens, readvar_tmp)
         155             :     !
         156             :     ! If field is on file:
         157             :     !
         158      146688 :     if (readvar_tmp) then
         159      144384 :         if (debug .and. masterproc) then
         160           0 :           write(iulog, '(2a,5(i0,a))') trim(subname),': field(',              &
         161           0 :                dim1b,':',dim1e,',',dim2b,':',dim2e, '), file(',dimlens(1),')'
         162           0 :           call shr_sys_flush(iulog)
         163             :         end if
         164             :       !
         165             :       ! Get array dimension id's and sizes
         166             :       !
         167      144384 :       arraydimsize(1) = (dim1e - dim1b + 1)
         168      144384 :       arraydimsize(2) = (dim2e - dim2b + 1)
         169      433152 :       do j = 1, 2
         170      433152 :         if (arraydimsize(j) /= size(field, j)) then
         171           0 :           write(errormsg, *) ': Mismatch between array bounds and field size for ', &
         172           0 :                trim(varname), ', dimension', j
         173           0 :           call endrun(trim(subname)//errormsg)
         174             :         end if
         175             :       end do
         176             :
         177      144384 :       if (ndims > 2) then
         178           0 :         call endrun(trim(subname)//': too many dimensions for '//trim(varname))
         179      144384 :       else if (ndims < 1) then
         180           0 :         call endrun(trim(subname)//': too few dimensions for '//trim(varname))
         181             :       else
         182             :          ! Check that the number of columns in the file matches the number of
         183             :          ! columns in the grid object.
         184      144384 :          if (dimlens(1) /= grid_dimlens(1) .and. .not. single_column) then
         185           0 :             readvar = .false.
         186           0 :             return
         187             :          end if
         188             :
         189             :         ! Check to make sure that the second dimension is time
         190      144384 :         if (ndims == 2) then
         191      135936 :           ierr = pio_inq_dimname(ncid, dimids(2), tmpname)
         192      135936 :           if (trim(tmpname) /= 'time') then
         193           0 :             call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
         194             :           end if
         195             :         end if
         196             :       end if
         197             :
         198      144384 :       if(ndims == 2) then
         199      135936 :         if(present(timelevel)) then
         200      135168 :           call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
         201             :         else
         202         768 :           call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
         203             :         end if
         204      135936 :         ndims = ndims - 1
         205             :       end if
         206             :
         207             :       ! nb: strt and cnt were initialized to 1
         208             :       ! all distributed array processing
         209           0 :       call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:ndims),    &
         210      144384 :            pio_double, iodesc)
         211      144384 :       call pio_read_darray(ncid, varid, iodesc, field, ierr)
         212      144384 :       if (present(fillvalue)) then
         213        5376 :          ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
         214             :       end if
         215             :
         216      144384 :       if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
         217             :
         218             :     end if  ! end of readvar_tmp
         219             :
         220      146688 :     readvar = readvar_tmp
         221             :
         222      146688 :     return
         223             :
         224      146688 :   end subroutine infld_real_1d_2d
         225             :
         226             :   !-----------------------------------------------------------------------
         227             :   !BOP
         228             :   !
         229             :   ! !IROUTINE: infld_real_2d_2d
         230             :   !
         231             :   ! !INTERFACE:
         232      141312 :   subroutine infld_real_2d_2d(varname, ncid, dimname1, dimname2,              &
         233      141312 :        dim1b, dim1e, dim2b, dim2e, field, readvar, gridname, timelevel,    &
         234             :        fillvalue)
         235             :     !
         236             :     ! !DESCRIPTION:
         237             :     ! Netcdf I/O of initial real field from netCDF file
         238             :     ! Read a 2-D field (or slice) into a 2-D variable
         239             :     !
         240             :     ! !USES
         241             :     !
         242             :
         243      146688 :     use pio,              only: pio_get_var, pio_read_darray, pio_setdebuglevel
         244             :     use pio,              only: PIO_MAX_NAME, pio_inq_dimname
         245             :     use cam_pio_utils,    only: cam_permute_array, calc_permutation
         246             :     use cam_pio_utils,    only: cam_pio_check_var, cam_pio_inq_var_fill
         247             :
         248             :     !
         249             :     ! !ARGUMENTS:
         250             :     implicit none
         251             :     character(len=*),  intent(in)     :: varname  ! variable name
         252             :     type(file_desc_t), intent(inout)  :: ncid     ! input unit
         253             :     character(len=*),  intent(in)     :: dimname1 ! name of 1st array dimensions of field on file (array order)
         254             :     character(len=*),  intent(in)     :: dimname2 ! name of 2nd array dimensions of field on file (array order)
         255             :     integer,           intent(in)     :: dim1b    ! start of first  dimension of array to be returned
         256             :     integer,           intent(in)     :: dim1e    ! end   of first  dimension of array to be returned
         257             :     integer,           intent(in)     :: dim2b    ! start of second dimension of array to be returned
         258             :     integer,           intent(in)     :: dim2e    ! end   of second dimension of array to be returned
         259             :     real(r8), target,  intent(out)    :: field(dim1b:dim1e,dim2b:dim2e) ! array to be returned (decomposed or global)
         260             :     logical,           intent(out)    :: readvar  ! true => variable is on initial dataset
         261             :     character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
         262             :     integer, optional, intent(in)     :: timelevel
         263             :     real(r8), optional, intent(out)   :: fillvalue
         264             :     !
         265             :     !EOP
         266             :     !
         267             :     ! !LOCAL VARIABLES:
         268             :     type(io_desc_t), pointer  :: iodesc
         269             :     integer                   :: grid_id   ! grid ID for data mapping
         270             :     integer                   :: i, j      ! indices
         271             :     integer                   :: ierr      ! error status
         272             :     type(var_desc_t)          :: varid     ! variable id
         273             :
         274             :     integer                   :: arraydimsize(2) ! field dimension lengths
         275             :     integer                   :: arraydimids(2) ! Dimension IDs
         276             :     integer                   :: permutation(2)
         277             :     logical                   :: ispermuted
         278             :
         279             :     integer                   :: ndims ! number of dimensions on file
         280             :     integer                   :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
         281             :     integer                   :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
         282             :
         283             :     ! Offsets for reading global variables
         284             :     integer                   :: strt(2) ! start lon, lat indices for netcdf 2-d
         285             :     integer                   :: cnt (2) ! lon, lat counts for netcdf 2-d
         286             :     character(len=PIO_MAX_NAME) :: tmpname
         287             :     character(len=128)        :: errormsg
         288             :
         289      141312 :     real(r8), pointer         :: tmp2d(:,:) ! input data for permutation
         290             :
         291             :     logical                   :: readvar_tmp ! if true, variable is on tape
         292             :     character(len=*), parameter :: subname='INFLD_REAL_2D_2D' ! subroutine name
         293             :     character(len=PIO_MAX_NAME) :: field_dnames(2)
         294             :     character(len=max_hcoordname_len) :: vargridname ! Name of variable's grid
         295             :
         296             :     ! For SCAM
         297             :     real(r8)                  :: closelat, closelon
         298             :     integer                   :: lonidx, latidx
         299             :
         300      141312 :     nullify(iodesc)
         301             :
         302             :     !
         303             :     !-----------------------------------------------------------------------
         304             :     !
         305             :     !    call pio_setdebuglevel(3)
         306             :
         307             :     ! Should we be using a different interface?
         308      141312 :     if ((trim(dimname1) == trim(dimname2)) .or. (len_trim(dimname2) == 0)) then
         309             :       call infld(varname, ncid, dimname1, dim1b, dim1e, dim2b, dim2e,         &
         310      141312 :            field, readvar, gridname, timelevel)
         311             :     else
         312             :
         313             :       !
         314             :       ! Error conditions
         315             :       !
         316           0 :       if (present(gridname)) then
         317           0 :         vargridname=trim(gridname)
         318             :       else
         319           0 :         vargridname='physgrid'
         320             :       end if
         321             :
         322           0 :       if (single_column .and. vargridname=='physgrid') then
         323           0 :          vargridname='physgrid_scm'
         324             :       end if
         325             :
         326           0 :       grid_id = cam_grid_id(trim(vargridname))
         327           0 :       if (.not. cam_grid_check(grid_id)) then
         328           0 :         if(masterproc) then
         329           0 :            write(errormsg, *)': invalid gridname, "',trim(vargridname),'", specified for field ',trim(varname)
         330             :         end if
         331           0 :         call endrun(trim(subname)//errormsg)
         332             :       end if
         333             :
         334           0 :       if (debug .and. masterproc) then
         335           0 :          write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(vargridname)
         336           0 :          call shr_sys_flush(iulog)
         337             :       end if
         338             :       !
         339             :       ! Read netCDF file
         340             :       !
         341             :       !
         342             :       ! Check if field is on file; get netCDF variable id
         343             :       !
         344           0 :       call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens, readvar_tmp)
         345             :       !
         346             :       ! If field is on file:
         347             :       !
         348           0 :       if (readvar_tmp) then
         349           0 :         if (debug .and. masterproc) then
         350           0 :           write(iulog, '(2a,6(i0,a))') trim(subname),': field(',              &
         351           0 :                dim1b,':',dim1e,',',dim2b,':',dim2e,                           &
         352           0 :                '), file(',dimlens(1),',',dimlens(2),')'
         353           0 :           call shr_sys_flush(iulog)
         354             :         end if
         355             :         !
         356             :         ! Get array dimension id's and sizes
         357             :         !
         358           0 :         ierr = PIO_inq_dimid(ncid, dimname1, arraydimids(1))
         359           0 :         ierr = PIO_inq_dimid(ncid, dimname2, arraydimids(2))
         360           0 :         arraydimsize(1) = (dim1e - dim1b + 1)
         361           0 :         arraydimsize(2) = (dim2e - dim2b + 1)
         362           0 :         do j = 1, 2
         363           0 :           if (arraydimsize(j) /= size(field, j)) then
         364           0 :             write(errormsg, *) ': Mismatch between array bounds and field size for ', &
         365           0 :                  trim(varname), ', dimension', j
         366           0 :             call endrun(trim(subname)//errormsg)
         367             :           end if
         368             :         end do
         369             :
         370           0 :         if (ndims > 3) then
         371           0 :           call endrun(trim(subname)//': too many dimensions for '//trim(varname))
         372           0 :         else if (ndims < 2) then
         373           0 :           call endrun(trim(subname)//': too few dimensions for '//trim(varname))
         374             :         else
         375             :           ! Check to make sure that the third dimension is time
         376           0 :           if (ndims == 3) then
         377           0 :             ierr = pio_inq_dimname(ncid, dimids(3), tmpname)
         378           0 :             if (trim(tmpname) /= 'time') then
         379           0 :               call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
         380             :             end if
         381             :           end if
         382             :         end if
         383             :
         384           0 :         if(ndims == 3) then
         385           0 :           if(present(timelevel)) then
         386           0 :             call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
         387             :           else
         388           0 :             call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
         389             :           end if
         390             :         end if
         391             :
         392           0 :         field_dnames(1) = dimname1
         393           0 :         field_dnames(2) = dimname2
         394           0 :         if (single_column) then
         395             :           ! This could be generalized but for now only handles a single point
         396           0 :           strt(1) = dim1b
         397           0 :           strt(2) = dim2b
         398           0 :           cnt = arraydimsize
         399           0 :           call shr_scam_getCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
         400           0 :           if (trim(field_dnames(1)) == 'lon') then
         401           0 :             strt(1) = lonidx ! First dim always lon for Eulerian dycore
         402             :           else
         403           0 :             call endrun(trim(subname)//': lon should be first dimension for '//trim(varname))
         404             :           end if
         405           0 :           if (trim(field_dnames(2)) == 'lat') then
         406           0 :             strt(2) = latidx
         407             :           else
         408           0 :             call endrun(trim(subname)//': lat dimension not found for '//trim(varname))
         409             :           end if
         410             :
         411             :           ! Check for permuted dimensions ('out of order' array)
         412           0 :           call calc_permutation(dimids, arraydimids, permutation, ispermuted)
         413           0 :           if (ispermuted) then
         414           0 :             call cam_permute_array(strt, permutation)
         415           0 :             call cam_permute_array(cnt, permutation)
         416           0 :             allocate(tmp2d(1:cnt(1), 1:cnt(2)))
         417           0 :             ierr = pio_get_var(ncid, varid, strt, cnt, tmp2d)
         418           0 :             do j = dim2b, dim2e
         419           0 :               do i = dim1b, dim1e
         420             :                 ! We don't need strt anymore, reuse it
         421           0 :                 strt(1) = i - dim1b + 1
         422           0 :                 strt(2) = j - dim2b + 1
         423           0 :                 call cam_permute_array(strt, permutation)
         424           0 :                 field(i,j) = tmp2d(strt(1), strt(2))
         425             :               end do
         426             :             end do
         427           0 :             deallocate(tmp2d)
         428             :           else
         429           0 :             ierr = pio_get_var(ncid, varid, strt, cnt, field)
         430             :           end if
         431             :         else
         432             :           ! All distributed array processing
         433             :           call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:2),      &
         434           0 :                pio_double, iodesc, field_dnames=field_dnames)
         435           0 :           call pio_read_darray(ncid, varid, iodesc, field, ierr)
         436           0 :           if (present(fillvalue)) then
         437           0 :              ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
         438             :           end if
         439             :         end if
         440             :
         441           0 :         if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
         442             :
         443             :       end if  ! end of readvar_tmp
         444             :
         445           0 :       readvar = readvar_tmp
         446             :
         447             :     end if ! end of call infld_real_1d_2d instead
         448             :
         449      282624 :   end subroutine infld_real_2d_2d
         450             :
         451             :
         452             :   !-----------------------------------------------------------------------
         453             :   !BOP
         454             :   !
         455             :   ! !IROUTINE: infld_real_2d_3d
         456             :   !
         457             :   ! !INTERFACE:
         458      125184 :   subroutine infld_real_2d_3d(varname, ncid, dimname1, dimname2,              &
         459             :        dim1b, dim1e, dim2b, dim2e, dim3b, dim3e,                              &
         460      125184 :        field, readvar, gridname, timelevel, fillvalue)
         461             :     !
         462             :     ! !DESCRIPTION:
         463             :     ! Netcdf I/O of initial real field from netCDF file
         464             :     ! Read a 2-D field (or slice) into a 3-D variable
         465             :     !
         466             :     ! !USES
         467             :     !
         468             :
         469      141312 :     use pio,              only: pio_get_var, pio_read_darray, pio_setdebuglevel
         470             :     use pio,              only: PIO_MAX_NAME, pio_inq_dimname
         471             :     use cam_pio_utils,    only: cam_pio_check_var, cam_pio_inq_var_fill
         472             :
         473             :     !
         474             :     ! !ARGUMENTS:
         475             :     implicit none
         476             :     character(len=*),  intent(in)     :: varname  ! variable name
         477             :     type(file_desc_t), intent(inout)  :: ncid     ! input unit
         478             :     character(len=*),  intent(in)     :: dimname1 ! name of 1st array dimensions of field on file (array order)
         479             :     character(len=*),  intent(in)     :: dimname2 ! name of 2nd array dimensions of field on file (array order)
         480             :     integer,           intent(in)     :: dim1b    ! start of first  dimension of array to be returned
         481             :     integer,           intent(in)     :: dim1e    ! end   of first  dimension of array to be returned
         482             :     integer,           intent(in)     :: dim2b    ! start of second dimension of array to be returned
         483             :     integer,           intent(in)     :: dim2e    ! end   of second dimension of array to be returned
         484             :     integer,           intent(in)     :: dim3b    ! start of third  dimension of array to be returned
         485             :     integer,           intent(in)     :: dim3e    ! end   of third  dimension of array to be returned
         486             :     real(r8), target,  intent(out)    :: field(dim1b:dim1e,dim2b:dim2e,dim3b:dim3e) ! array to be returned (decomposed or global)
         487             :     logical,           intent(out)    :: readvar  ! true => variable is on initial dataset
         488             :     character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
         489             :     integer, optional, intent(in)     :: timelevel
         490             :     real(r8), optional, intent(out)   :: fillvalue
         491             :     !
         492             :     !EOP
         493             :     !
         494             :     ! !LOCAL VARIABLES:
         495             :     type(io_desc_t), pointer  :: iodesc
         496             :     integer                   :: grid_id   ! grid ID for data mapping
         497             :     integer                   :: j         ! index
         498             :     integer                   :: ierr      ! error status
         499             :     type(var_desc_t)          :: varid     ! variable id
         500             :
         501             :     integer                   :: arraydimsize(3) ! field dimension lengths
         502             :
         503             :     integer                   :: ndims ! number of dimensions
         504             :     integer                   :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
         505             :     integer                   :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
         506             :     integer                   :: grid_dimlens(2)
         507             :
         508             :     ! Offsets for reading global variables
         509             :     integer                   :: strt(3) = 1 ! start ncol, lev indices for netcdf 2-d
         510             :     integer                   :: cnt (3) = 1 ! ncol, lev counts for netcdf 2-d
         511             :     character(len=PIO_MAX_NAME) :: tmpname
         512             :
         513             :     logical                   :: readvar_tmp ! if true, variable is on tape
         514             :     character(len=*), parameter :: subname='INFLD_REAL_2D_3D' ! subroutine name
         515             :     character(len=128)        :: errormsg
         516             :     character(len=PIO_MAX_NAME) :: field_dnames(2)
         517             :     character(len=PIO_MAX_NAME) :: file_dnames(3)
         518             :     character(len=max_hcoordname_len) :: vargridname ! Name of variable's grid
         519             :
         520             :     !
         521             :     !-----------------------------------------------------------------------
         522             :     !
         523             :     !    call pio_setdebuglevel(3)
         524             :
         525      125184 :     nullify(iodesc)
         526             :
         527             :     !
         528             :     ! Error conditions
         529             :     !
         530      125184 :     if (present(gridname)) then
         531      125184 :       vargridname=trim(gridname)
         532             :     else
         533           0 :       vargridname='physgrid'
         534             :     end if
         535             :
         536             :     ! if running single column mode then we need to use scm grid to read proper column
         537      125184 :     if (single_column .and. vargridname=='physgrid') then
         538           0 :        vargridname='physgrid_scm'
         539             :     end if
         540             :
         541      125184 :     grid_id = cam_grid_id(trim(vargridname))
         542      125184 :     if (.not. cam_grid_check(grid_id)) then
         543           0 :       if(masterproc) then
         544           0 :          write(errormsg, *)': invalid gridname, "',trim(vargridname),'", specified for field ',trim(varname)
         545             :       end if
         546           0 :       call endrun(trim(subname)//errormsg)
         547             :     end if
         548             :
         549      125184 :     if (debug .and. masterproc) then
         550           0 :       write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(vargridname)
         551           0 :       call shr_sys_flush(iulog)
         552             :     end if
         553             :
         554             :     ! Get the number of columns in the global grid.
         555      125184 :     call cam_grid_dimensions(grid_id, grid_dimlens)
         556             :     !
         557             :     ! Read netCDF file
         558             :     !
         559             :     !
         560             :     ! Check if field is on file; get netCDF variable id
         561             :     !
         562             :     call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens, &
         563      125184 :        readvar_tmp, dimnames=file_dnames)
         564             :
         565             :     ! If field is on file:
         566             :     !
         567      125184 :     if (readvar_tmp) then
         568      121344 :       if (debug .and. masterproc) then
         569           0 :         write(iulog, '(2a,8(i0,a))') trim(subname),': field(',                &
         570           0 :              dim1b,':',dim1e,',',dim2b,':',dim2e,',',dim3b,':',dim3e,         &
         571           0 :              '), file(',dimlens(1),',',dimlens(2),')'
         572           0 :         call shr_sys_flush(iulog)
         573             :       end if
         574             :       !
         575             :       ! Get array dimension id's and sizes
         576             :       !
         577      121344 :       arraydimsize(1) = (dim1e - dim1b + 1)
         578      121344 :       arraydimsize(2) = (dim2e - dim2b + 1)
         579      121344 :       arraydimsize(3) = (dim3e - dim3b + 1)
         580      485376 :       do j = 1, 3
         581      485376 :         if (arraydimsize(j) /= size(field, j)) then
         582           0 :           write(errormsg, *) ': Mismatch between array bounds and field size for ', &
         583           0 :                trim(varname), ', dimension', j
         584           0 :           call endrun(trim(subname)//errormsg)
         585             :         end if
         586             :       end do
         587             :
         588      121344 :       if (ndims > 3) then
         589           0 :         call endrun(trim(subname)//': too many dimensions for '//trim(varname))
         590      121344 :       else if (ndims < 2) then
         591           0 :         call endrun(trim(subname)//': too few dimensions for '//trim(varname))
         592             :       else
         593             :          ! Check that the number of columns in the file matches the number of
         594             :          ! columns in the grid object.
         595      121344 :          if (dimlens(1) /= grid_dimlens(1) .and. dimlens(2) /= grid_dimlens(1) .and. .not. single_column) then
         596           0 :             readvar = .false.
         597           0 :             return
         598             :          end if
         599             :
         600             :         ! Check to make sure that the 3rd dimension is time
         601      121344 :         if (ndims == 3) then
         602      112128 :           ierr = pio_inq_dimname(ncid, dimids(3), tmpname)
         603      224256 :           if (to_lower(trim(tmpname)) /= 'time') then
         604      112128 :             call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
         605             :           end if
         606             :         end if
         607             :       end if
         608             :
         609      121344 :       if(ndims == 3) then
         610      112128 :         if(present(timelevel)) then
         611       78336 :           call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
         612             :         else
         613       33792 :           call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
         614             :         end if
         615      112128 :         ndims = ndims - 1
         616             :       end if
         617             :
         618      121344 :       field_dnames(1) = dimname1
         619      121344 :       field_dnames(2) = dimname2
         620             :       ! NB: strt and cnt were initialized to 1
         621             :       ! All distributed array processing
         622             :       call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:2),         &
         623             :            pio_double, iodesc, field_dnames=field_dnames,                   &
         624      121344 :            file_dnames=file_dnames(1:2))
         625      121344 :       call pio_read_darray(ncid, varid, iodesc, field, ierr)
         626      121344 :       if (present(fillvalue)) then
         627       33792 :          ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
         628             :       end if
         629             :
         630      121344 :       if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
         631             :
         632             :     end if  ! end of readvar_tmp
         633             :
         634      125184 :     readvar = readvar_tmp
         635             :
         636      125184 :     return
         637             :
         638      250368 :   end subroutine infld_real_2d_3d
         639             :
         640             :   !-----------------------------------------------------------------------
         641             :   !BOP
         642             :   !
         643             :   ! !IROUTINE: infld_real_3d_3d
         644             :   !
         645             :   ! !INTERFACE:
         646       57600 :   subroutine infld_real_3d_3d(varname, ncid, dimname1, dimname2, dimname3,    &
         647             :        dim1b, dim1e, dim2b, dim2e, dim3b, dim3e,                              &
         648       57600 :        field, readvar, gridname, timelevel, fillvalue)
         649             :     !
         650             :     ! !DESCRIPTION:
         651             :     ! Netcdf I/O of initial real field from netCDF file
         652             :     ! Read a 3-D field (or slice) into a 3-D variable
         653             :     !
         654             :     ! !USES
         655             :     !
         656             :
         657      125184 :     use pio,              only: pio_get_var, pio_read_darray, pio_setdebuglevel
         658             :     use pio,              only: PIO_MAX_NAME, pio_inq_dimname
         659             :     use cam_pio_utils,    only: cam_permute_array, calc_permutation
         660             :     use cam_pio_utils,    only: cam_pio_check_var, cam_pio_inq_var_fill
         661             :
         662             :     !
         663             :     ! !ARGUMENTS:
         664             :     implicit none
         665             :     character(len=*),  intent(in)     :: varname  ! variable name
         666             :     type(file_desc_t), intent(inout)  :: ncid     ! input unit
         667             :     character(len=*),  intent(in)     :: dimname1 ! name of 1st array dimensions of field on file (array order)
         668             :     character(len=*),  intent(in)     :: dimname2 ! name of 2nd array dimensions of field on file (array order)
         669             :     character(len=*),  intent(in)     :: dimname3 ! name of 3rd array dimensions of field on file (array order)
         670             :     integer,           intent(in)     :: dim1b    ! start of first  dimension of array to be returned
         671             :     integer,           intent(in)     :: dim1e    ! end   of first  dimension of array to be returned
         672             :     integer,           intent(in)     :: dim2b    ! start of second dimension of array to be returned
         673             :     integer,           intent(in)     :: dim2e    ! end   of second dimension of array to be returned
         674             :     integer,           intent(in)     :: dim3b    ! start of third  dimension of array to be returned
         675             :     integer,           intent(in)     :: dim3e    ! end   of third  dimension of array to be returned
         676             :     real(r8), target,  intent(out)    :: field(dim1b:dim1e,dim2b:dim2e,dim3b:dim3e) ! array to be returned (decomposed or global)
         677             :     logical,           intent(out)    :: readvar  ! true => variable is on initial dataset
         678             :     character(len=*), optional, intent(in) :: gridname ! Name of variable's grid
         679             :     integer, optional, intent(in)     :: timelevel
         680             :     real(r8), optional, intent(out)   :: fillvalue
         681             :     !
         682             :     !EOP
         683             :     !
         684             :     ! !LOCAL VARIABLES:
         685             :     type(io_desc_t), pointer  :: iodesc
         686             :     integer                   :: grid_id   ! grid ID for data mapping
         687             :     integer                   :: i, j, k   ! indices
         688             :     integer                   :: ierr      ! error status
         689             :     type(var_desc_t)          :: varid     ! variable id
         690             :
         691             :     integer                   :: arraydimsize(3) ! field dimension lengths
         692             :     integer                   :: arraydimids(3) ! Dimension IDs
         693             :     integer                   :: permutation(3)
         694             :     logical                   :: ispermuted
         695             :
         696             :     integer                   :: ndims ! number of dimensions
         697             :     integer                   :: pdims ! number of dimensions w/o timeslice
         698             :     integer                   :: dimids(PIO_MAX_VAR_DIMS) ! file variable dims
         699             :     integer                   :: dimlens(PIO_MAX_VAR_DIMS) ! file variable shape
         700             :
         701             :     ! Offsets for reading global variables
         702             :     integer                   :: strt(3)     ! start lon, lev, lat indices for netcdf 3-d
         703             :     integer                   :: cnt (3)     ! lon, lat counts for netcdf 3-d
         704             :     character(len=PIO_MAX_NAME) :: tmpname
         705             :
         706       57600 :     real(r8), pointer         :: tmp3d(:,:,:) ! input data for permutation
         707             :
         708             :     logical                   :: readvar_tmp ! if true, variable is on tape
         709             :     character(len=*), parameter :: subname='INFLD_REAL_3D_3D' ! subroutine name
         710             :     character(len=128)        :: errormsg
         711             :     character(len=PIO_MAX_NAME) :: field_dnames(3)
         712             :     character(len=PIO_MAX_NAME) :: file_dnames(4)
         713             :     character(len=max_hcoordname_len) :: vargridname ! Name of variable's grid
         714             :
         715             :     ! For SCAM
         716             :     real(r8)                  :: closelat, closelon
         717             :     integer                   :: lonidx, latidx
         718             :
         719       57600 :     nullify(iodesc)
         720             :
         721             :     !
         722             :     !-----------------------------------------------------------------------
         723             :     !
         724             :     !    call pio_setdebuglevel(3)
         725             :
         726             :     ! Should we be using a different interface?
         727       57600 :     if ((trim(dimname1) == trim(dimname2)) .or. (len_trim(dimname2) == 0)) then
         728             :       call infld(varname, ncid, dimname1, dimname3,                           &
         729             :            dim1b, dim1e, dim2b, dim2e, dim3b, dim3e,                          &
         730           0 :            field, readvar, gridname, timelevel)
         731       57600 :     else if ((trim(dimname1) == trim(dimname3)) .or. (len_trim(dimname3) == 0)) then
         732             :       call infld(varname, ncid, dimname1, dimname2,                           &
         733             :            dim1b, dim1e, dim2b, dim2e, dim3b, dim3e,                          &
         734       57600 :            field, readvar, gridname, timelevel)
         735             :     else
         736             :       !
         737             :       ! Error conditions
         738             :       !
         739           0 :       if (present(gridname)) then
         740           0 :         vargridname=trim(gridname)
         741             :       else
         742           0 :         vargridname='physgrid'
         743             :       end if
         744             :
         745             :       ! if running single column mode then we need to use scm grid to read proper column
         746           0 :       if (single_column .and. vargridname=='physgrid') then
         747           0 :          vargridname='physgrid_scm'
         748             :       end if
         749             :
         750           0 :       grid_id = cam_grid_id(trim(vargridname))
         751           0 :       if (.not. cam_grid_check(grid_id)) then
         752           0 :         if(masterproc) then
         753           0 :            write(errormsg, *)': invalid gridname, "',trim(vargridname),'", specified for field ',trim(varname)
         754             :         end if
         755           0 :         call endrun(trim(subname)//errormsg)
         756             :       end if
         757             :
         758           0 :       if (debug .and. masterproc) then
         759           0 :         write(iulog, '(5a)') trim(subname),': field = ',trim(varname),', grid = ',trim(vargridname)
         760           0 :         call shr_sys_flush(iulog)
         761             :       end if
         762             :       !
         763             :       ! Read netCDF file
         764             :       !
         765             :       !
         766             :       ! Check if field is on file; get netCDF variable id
         767             :       !
         768             :       call cam_pio_check_var(ncid, varname, varid, ndims, dimids, dimlens,    &
         769           0 :            readvar_tmp, dimnames=file_dnames)
         770             :       !
         771             :       ! If field is on file:
         772             :       !
         773           0 :       if (readvar_tmp) then
         774           0 :         if (debug .and. masterproc) then
         775           0 :           write(iulog, '(2a,9(i0,a))') trim(subname),': field(',              &
         776           0 :                dim1b,':',dim1e,',',dim2b,':',dim2e,',',dim3b,':',dim3e,       &
         777           0 :                '), file(',dimlens(1),',',dimlens(2),',',dimlens(3),')'
         778           0 :           call shr_sys_flush(iulog)
         779             :         end if
         780             :         !
         781             :         ! Get array dimension id's and sizes
         782             :         !
         783           0 :         ierr = PIO_inq_dimid(ncid, dimname1, arraydimids(1))
         784           0 :         ierr = PIO_inq_dimid(ncid, dimname2, arraydimids(2))
         785           0 :         ierr = PIO_inq_dimid(ncid, dimname3, arraydimids(3))
         786           0 :         arraydimsize(1) = (dim1e - dim1b + 1)
         787           0 :         arraydimsize(2) = (dim2e - dim2b + 1)
         788           0 :         arraydimsize(3) = (dim3e - dim3b + 1)
         789             :
         790           0 :         do j = 1, 3
         791           0 :           if (arraydimsize(j) /= size(field, j)) then
         792           0 :             write(errormsg, *) ': Mismatch between array bounds and field size for ', &
         793           0 :                  trim(varname), ', dimension', j
         794           0 :             call endrun(trim(subname)//errormsg)
         795             :           end if
         796             :         end do
         797             :
         798           0 :         pdims = ndims
         799           0 :         if (ndims > 4) then
         800           0 :           call endrun(trim(subname)//': too many dimensions for '//trim(varname))
         801           0 :         else if (ndims < 3) then
         802           0 :           call endrun(trim(subname)//': too few dimensions for '//trim(varname))
         803             :         else
         804             :           ! Check to make sure that the fourth dimension is time
         805           0 :           if (ndims == 4) then
         806           0 :             ierr = pio_inq_dimname(ncid, dimids(4), tmpname)
         807           0 :             if (trim(tmpname) /= 'time') then
         808           0 :               call endrun(trim(subname)//': dimension mismatch for '//trim(varname))
         809             :             end if
         810           0 :             pdims = 3
         811             :           end if
         812             :         end if
         813             :
         814           0 :         if(ndims == 4) then
         815           0 :           if(present(timelevel)) then
         816           0 :             call pio_setframe(ncid, varid, int(timelevel,kind=pio_offset_kind))
         817             :           else
         818           0 :             call pio_setframe(ncid, varid, int(1,kind=pio_offset_kind))
         819             :           end if
         820             :         end if
         821             :
         822           0 :         field_dnames(1) = dimname1
         823           0 :         field_dnames(2) = dimname2
         824           0 :         field_dnames(3) = dimname3
         825             :
         826           0 :         if (single_column) then
         827             :           ! This could be generalized but for now only handles a single point
         828           0 :           strt(1) = dim1b
         829           0 :           strt(2) = dim2b
         830           0 :           strt(3) = dim3b
         831           0 :           cnt = arraydimsize
         832           0 :           call shr_scam_getCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
         833           0 :           if (trim(field_dnames(1)) == 'lon') then
         834           0 :             strt(1) = lonidx ! First dim always lon for Eulerian dycore
         835             :           else
         836           0 :             call endrun(trim(subname)//': lon should be first dimension for '//trim(varname))
         837             :           end if
         838           0 :           if (trim(field_dnames(2)) == 'lat') then
         839           0 :             strt(2) = latidx
         840           0 :           else if (trim(field_dnames(3)) == 'lat') then
         841           0 :             strt(3) = latidx
         842             :           else
         843           0 :             call endrun(trim(subname)//': lat dimension not found for '//trim(varname))
         844             :           end if
         845             :
         846             :           ! Check for permuted dimensions ('out of order' array)
         847           0 :           call calc_permutation(dimids, arraydimids, permutation, ispermuted)
         848           0 :           if (ispermuted) then
         849           0 :             call cam_permute_array(strt, permutation)
         850           0 :             call cam_permute_array(cnt, permutation)
         851           0 :             allocate(tmp3d(1:cnt(1), 1:cnt(2), 1:cnt(3)))
         852           0 :             ierr = pio_get_var(ncid, varid, strt, cnt, tmp3d)
         853           0 :             do k = dim3b, dim3e
         854           0 :               do j = dim2b, dim2e
         855           0 :                 do i = dim1b, dim1e
         856             :                   ! We don't need strt anymore, reuse it
         857           0 :                   strt(1) = i - dim1b + 1
         858           0 :                   strt(2) = j - dim2b + 1
         859           0 :                   strt(3) = k - dim3b + 1
         860           0 :                   call cam_permute_array(strt, permutation)
         861           0 :                   field(i,j,k) = tmp3d(strt(1), strt(2), strt(3))
         862             :                 end do
         863             :               end do
         864             :             end do
         865           0 :             deallocate(tmp3d)
         866             :           else
         867           0 :             ierr = pio_get_var(ncid, varid, strt, cnt, field)
         868             :           end if
         869             :         else
         870             :           ! All distributed array processing
         871           0 :           call cam_grid_get_decomp(grid_id, arraydimsize, dimlens(1:pdims),   &
         872             :                pio_double, iodesc, field_dnames=field_dnames,                 &
         873           0 :                file_dnames=file_dnames(1:3))
         874           0 :           call pio_read_darray(ncid, varid, iodesc, field, ierr)
         875           0 :           if (present(fillvalue)) then
         876           0 :              ierr = cam_pio_inq_var_fill(ncid, varid, fillvalue)
         877             :           end if
         878             :         end if ! end of single column
         879             :
         880           0 :         if (masterproc) write(iulog,*) subname//': read field '//trim(varname)
         881             :
         882             :       end if  ! end of readvar_tmp
         883             :
         884           0 :       readvar = readvar_tmp
         885             :
         886             :     end if ! end of call infld_real_2d_3d instead
         887             :
         888      115200 :   end subroutine infld_real_3d_3d
         889             :
         890             :
         891             : end module ncdio_atm

Generated by: LCOV version 1.14