LCOV - code coverage report
Current view: | top level - control - cam_history_support.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 221 | 705 | 31.3 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 23 | 55 | 41.8 %

              Line data    Source code
           1             : module cam_history_support
           2             :
           3             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           4             : !!
           5             : !!  cam_history_support is used by cam_history as well as by the dycores
           6             : !!    (for vertical coordinate and "mdim" support). Some parameters are
           7             : !!    also referenced by cam_grid_support (although those could be copied
           8             : !!    if necessary).
           9             : !!
          10             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          11             :
          12             :   use shr_kind_mod,     only: r8=>shr_kind_r8, shr_kind_cl, shr_kind_cxx
          13             :   use pio,              only: var_desc_t, file_desc_t
          14             :   use cam_abortutils,   only: endrun
          15             :   use cam_logfile,      only: iulog
          16             :   use spmd_utils,       only: masterproc
          17             :   use cam_grid_support, only: cam_grid_patch_t, cam_grid_header_info_t
          18             :   use cam_grid_support, only: max_hcoordname_len, maxsplitfiles
          19             :   use cam_pio_utils,    only: cam_pio_handle_error
          20             :
          21             :   implicit none
          22             :   private
          23             :   save
          24             :
          25             :   integer, parameter, public :: max_string_len = shr_kind_cxx
          26             :   integer, parameter, public :: max_chars = shr_kind_cl         ! max chars for char variables
          27             :   integer, parameter, public :: field_op_len  = 3               ! max chars for field operation string (sum/dif)
          28             :   integer, parameter, public :: fieldname_len = 32              ! max chars for field name
          29             :   integer, parameter, public :: fieldname_suffix_len =  3       ! length of field name suffix ("&IC")
          30             :   integer, parameter, public :: fieldname_lenp2 = fieldname_len + 2 ! allow for extra characters
          31             :   ! max_fieldname_len = max chars for field name (including suffix)
          32             :   integer, parameter, public :: max_fieldname_len    = fieldname_len + fieldname_suffix_len
          33             :
          34             :   integer, parameter, public :: pflds  = 1000      ! max number of fields for namelist entries fincl and fexcl
          35             :                                                    ! also used in write restart
          36             :   integer, parameter, public :: ptapes = 12        ! max number of tapes
          37             :
          38             :   ! A special symbol for declaring a field which has no vertical or
          39             :   ! non-grid dimensions. It is here (rather than cam_history) so that it
          40             :   ! be checked by add_hist_coord
          41             :   character(len=10), parameter, public :: horiz_only = 'horiz_only'
          42             :
          43             :   type, public :: history_patch_t
          44             :     character(len=max_chars)          :: namelist_entry =  ''
          45             :     ! lon_axis_name and lat_axis_name are not used if collected_output = .true.
          46             :     character(len=max_fieldname_len)  :: lon_axis_name       =  ''
          47             :     character(len=max_fieldname_len)  :: lat_axis_name       =  ''
          48             :     logical                           :: collected_output
          49             :     ! There is one patch for every grid and one header_info for every patch
          50             :     type(cam_grid_patch_t),        pointer :: patches(:)     => NULL()
          51             :     type (cam_grid_header_info_t), pointer :: header_info(:) => NULL()
          52             :   contains
          53             :     procedure                   :: write_attrs  => history_patch_write_attrs
          54             :     procedure                   :: write_vals   => history_patch_write_vals
          55             :     procedure                   :: field_name   => history_patch_field_name
          56             :     procedure                   :: num_hdims    => history_patch_num_hdims
          57             :     procedure                   :: get_var_data => history_patch_get_var_data
          58             :     procedure                   :: write_var    => history_patch_write_var
          59             :     procedure                   :: compact      => history_patch_compact
          60             :     procedure                   :: active_cols  => history_patch_active_cols
          61             :     procedure                   :: deallocate   => history_patch_deallocate
          62             :   end type history_patch_t
          63             :
          64             : !
          65             : ! dim_index_2d, dim_index_3d: 2-D & 3-D dimension index lower & upper bounds
          66             : !
          67             :   type, public :: dim_index_2d        ! 2-D dimension index
          68             :      integer :: beg1, end1            ! lower & upper bounds of 1st dimension
          69             :      integer :: beg2, end2            ! lower & upper bounds of 2nd dimension
          70             :    contains
          71             :      procedure :: dim_sizes_2d  => dim_index_2d_dim_sizes_2d
          72             :      procedure :: dim_sizes_arr => dim_index_2d_dim_size_arr
          73             :      generic   :: dim_sizes     => dim_sizes_arr, dim_sizes_2d
          74             :   end type dim_index_2d
          75             :
          76             :   type, public :: dim_index_3d       ! 3-D dimension index
          77             :     integer :: beg1, end1            ! lower & upper bounds of 1st dimension
          78             :     integer :: beg2, end2            ! lower & upper bounds of 2nd dimension
          79             :     integer :: beg3, end3            ! lower & upper bounds of 3rd dimension
          80             :    contains
          81             :      procedure :: dim_sizes_3d  => dim_index_3d_dim_sizes_3d
          82             :      procedure :: dim_sizes_arr => dim_index_3d_dim_size_arr
          83             :      generic   :: dim_sizes     => dim_sizes_arr, dim_sizes_3d
          84             :   end type dim_index_3d
          85             :
          86             :   !---------------------------------------------------------------------------
          87             :   !
          88             :   ! field_info: A derived type containing information in an addfld call.
          89             :   !
          90             :   !---------------------------------------------------------------------------
          91             :   type, public :: field_info
          92             :
          93             :     logical :: flag_xyfill                   ! non-applicable xy points flagged with fillvalue
          94             :     logical :: is_subcol                     ! .true. iff field output as subcol
          95             :     integer, pointer :: mdims(:) => NULL()   ! indicies into hist_coords list
          96             :     integer, pointer :: shape(:) => NULL()   ! shape of field on file
          97             :
          98             :     real(r8) :: fillvalue                    ! fillvalue for this variable, set to default if not explicit in addfld
          99             :
         100             :     integer :: numlev                        ! vertical dimension (.nc file and internal arr)
         101             :
         102             :     integer :: begdim1                       ! on-node dim1 start index
         103             :     integer :: enddim1                       ! on-node dim1 end index
         104             :
         105             :     integer :: begdim2                       ! on-node dim2 start index
         106             :     integer :: enddim2                       ! on-node dim2 end index
         107             :
         108             :     integer :: begdim3                       ! on-node chunk or lat start index
         109             :     integer :: enddim3                       ! on-node chunk or lat end index
         110             :
         111             :     logical :: colperchunk                   ! .true. iff ncols /= chunksize
         112             :
         113             :     integer :: decomp_type                   ! type of decomposition (e.g., physics or dynamics)
         114             :
         115             :     ! meridional and zonal complements are for fields defined as part of a
         116             :     ! 2-D vector. These IDs are used to facilitate interpolated history output
         117             :     ! At most one of these will be a positive field ID.
         118             :     integer :: meridional_complement         ! meridional field id or -1
         119             :     integer :: zonal_complement              ! zonal field id or -1
         120             :
         121             :     character(len=field_op_len) :: field_op = ''        ! 'sum' or 'dif'
         122             :     integer                     :: op_field1_id         ! first field id or -1
         123             :     integer                     :: op_field2_id         ! second field id or -1
         124             :
         125             :     character(len=max_fieldname_len) :: name ! field name
         126             :     character(len=max_chars) :: long_name    ! long name
         127             :     character(len=max_chars) :: units        ! units
         128             :     character(len=3)         :: mixing_ratio ! 'dry' or 'wet'
         129             :     character(len=max_chars) :: sampling_seq ! sampling sequence - if not every timestep, how often field is sampled
         130             :     ! (i.e., how often "outfld" is called):  every other; only during LW/SW
         131             :     ! radiation calcs; etc.
         132             :     character(len=max_chars) :: cell_methods ! optional cell_methods attribute
         133             :   contains
         134             :     procedure :: is_composed => field_info_is_composed
         135             :     procedure :: get_shape   => field_info_get_shape
         136             :     procedure :: get_bounds  => field_info_get_bounds
         137             :     procedure :: get_dims_2d => field_info_get_dims_2d
         138             :     procedure :: get_dims_3d => field_info_get_dims_3d
         139             :     generic   :: get_dims    => get_dims_2d, get_dims_3d
         140             :   end type field_info
         141             :
         142             :   real(r8), parameter, public :: fillvalue = 1.e36_r8     ! fill value for netcdf fields
         143             :
         144             :
         145             :   !---------------------------------------------------------------------------
         146             :   !
         147             :   ! hentry: elements of an entry in the list of active fields on a single
         148             :   !         history file
         149             :   !         nacs is an accumulation counter which normally counts an entire
         150             :   !         chunk (physics) or block (dynamics) as accumulated as a single
         151             :   !         entity. The per-chunk counting avoids counting multiple outfld
         152             :   !         calls as multiple accumulations. Only the value of the first chunk
         153             :   !         or block is written to or read from a history restart file.
         154             :   !         For certain actions (e.g., only accumulating on
         155             :   !         non-fillvalue or accumulating based on local time), nacs has an
         156             :   !         entry for every column.
         157             :   !         nacs does not keep track of levels
         158             :   !
         159             :   !---------------------------------------------------------------------------
         160             :   type, public:: hentry
         161             :     type (field_info)         :: field                 ! field information
         162             :     character(len=1)          :: avgflag               ! averaging flag
         163             :     character(len=max_chars)  :: time_op               ! time operator (e.g. max, min, avg)
         164             :     character(len=max_fieldname_len)  :: op_field1     ! field1 name for sum or dif operation
         165             :     character(len=max_fieldname_len)  :: op_field2     ! field2 name for sum or dif operation
         166             :
         167             :     integer                   :: hwrt_prec             ! history output precision
         168             :     real(r8),         pointer :: hbuf(:,:,:) => NULL()
         169             :     real(r8),         private :: hbuf_integral         ! area weighted integral of active field
         170             :     real(r8),         pointer :: sbuf(:,:,:) => NULL() ! for standard deviation
         171             :     real(r8),         pointer :: wbuf(:,:)   => NULL() ! pointer to area weights
         172             :     type(var_desc_t), pointer :: varid(:)    => NULL() ! variable ids
         173             :     integer,          pointer :: nacs(:,:)   => NULL() ! accumulation counter
         174             :     type(var_desc_t), pointer :: nacs_varid  => NULL()
         175             :     integer                   :: beg_nstep             ! starting time step for nstep normalization
         176             :     type(var_desc_t), pointer :: beg_nstep_varid=> NULL()
         177             :     type(var_desc_t), pointer :: sbuf_varid  => NULL()
         178             :     type(var_desc_t), pointer :: wbuf_varid  => NULL()
         179             :   contains
         180             :     procedure :: get_global   => hentry_get_global
         181             :     procedure :: put_global   => hentry_put_global
         182             :   end type hentry
         183             :
         184             :   !---------------------------------------------------------------------------
         185             :   !
         186             :   !  active_entry: derived type containing information for a history tape
         187             :   !
         188             :   !---------------------------------------------------------------------------
         189             :   type, public:: active_entry
         190             :
         191             :     type(hentry), pointer :: hlist(:)
         192             :
         193             :     integer,               pointer :: grid_ids(:) => NULL()
         194             :     type(history_patch_t), pointer :: patches(:)  => NULL()
         195             :
         196             :     !
         197             :     ! PIO ids
         198             :     !
         199             :
         200             :     type(file_desc_t) :: Files(maxsplitfiles)  ! PIO file ids
         201             :
         202             :     type(var_desc_t) :: mdtid            ! var id for timestep
         203             :     type(var_desc_t) :: ndbaseid         ! var id for base day
         204             :     type(var_desc_t) :: nsbaseid         ! var id for base seconds of base day
         205             :     type(var_desc_t) :: nbdateid         ! var id for base date
         206             :     type(var_desc_t) :: nbsecid          ! var id for base seconds of base date
         207             :     type(var_desc_t) :: ndcurid          ! var id for current day
         208             :     type(var_desc_t) :: nscurid          ! var id for current seconds of current day
         209             :     type(var_desc_t) :: dateid           ! var id for current date
         210             :     type(var_desc_t) :: co2vmrid         ! var id for co2 volume mixing ratio
         211             :     type(var_desc_t) :: ch4vmrid         ! var id for ch4 volume mixing ratio
         212             :     type(var_desc_t) :: n2ovmrid         ! var id for n2o volume mixing ratio
         213             :     type(var_desc_t) :: f11vmrid         ! var id for f11 volume mixing ratio
         214             :     type(var_desc_t) :: f12vmrid         ! var id for f12 volume mixing ratio
         215             :     type(var_desc_t) :: sol_tsiid        ! var id for total solar irradiance (W/m2)
         216             :     type(var_desc_t) :: datesecid        ! var id for curent seconds of current date
         217             : #if ( defined BFB_CAM_SCAM_IOP )
         218             :     type(var_desc_t) :: bdateid         ! var id for base date
         219             :     type(var_desc_t) :: tsecid        ! var id for curent seconds of current date
         220             : #endif
         221             :     type(var_desc_t) :: nstephid         ! var id for current timestep
         222             :     type(var_desc_t) :: timeid           ! var id for time
         223             :     type(var_desc_t) :: tbndid           ! var id for time_bounds
         224             :     type(var_desc_t) :: date_writtenid   ! var id for date time sample written
         225             :     type(var_desc_t) :: time_writtenid   ! var id for time time sample written
         226             :     type(var_desc_t) :: f107id           ! var id for f107
         227             :     type(var_desc_t) :: f107aid          ! var id for f107a
         228             :     type(var_desc_t) :: f107pid          ! var id for f107p
         229             :     type(var_desc_t) :: kpid             ! var id for kp
         230             :     type(var_desc_t) :: apid             ! var id for ap
         231             :     type(var_desc_t) :: byimfid          ! var id IMF BY
         232             :     type(var_desc_t) :: bzimfid          ! var id IMF BZ
         233             :     type(var_desc_t) :: swvelid          ! var id solar wind velocity
         234             :     type(var_desc_t) :: swdenid          ! var id solar wind density
         235             :     type(var_desc_t) :: colat_crit1_id   ! critical colatitude
         236             :     type(var_desc_t) :: colat_crit2_id   ! critical colatitude
         237             :
         238             :   end type active_entry
         239             :
         240             :   !---------------------------------------------------------------------------
         241             :   !
         242             :   !  formula_terms_t: Information for formula terms (CF convention) variables
         243             :   !                   Used to add a formula-terms variable to the history file
         244             :   !                   Also adds a string, '<name>: <var_name>' to the parent
         245             :   !                   mdim's 'formula_terms' attribute.
         246             :   !
         247             :   !---------------------------------------------------------------------------
         248             :   type, public :: formula_terms_t
         249             :     character(len=max_fieldname_len) :: a_name = ''   ! 'A' term variable name
         250             :     character(len=max_string_len)    :: a_long_name = '' ! 'A' long name
         251             :     real(r8), pointer                :: a_values(:) => null() ! 'A' variable values
         252             :     character(len=max_fieldname_len) :: b_name = ''   ! 'B' term variable name
         253             :     character(len=max_string_len)    :: b_long_name = '' ! 'B' long name
         254             :     real(r8), pointer                :: b_values(:) => null() ! 'B' variable values
         255             :     character(len=max_fieldname_len) :: p0_name = ''  ! 'p0' term variable name
         256             :     character(len=max_string_len)    :: p0_long_name = '' ! 'p0' long name
         257             :     character(len=max_chars)         :: p0_units = '' ! 'p0' variable units
         258             :     real(r8)                         :: p0_value = fillvalue ! 'p0' variable values
         259             :     character(len=max_fieldname_len) :: ps_name = ''  ! 'ps' term variable name
         260             :   end type formula_terms_t
         261             :
         262             :   !---------------------------------------------------------------------------
         263             :   !
         264             :   !  hist_coord_t: Information for history variable dimension attributes
         265             :   !
         266             :   !---------------------------------------------------------------------------
         267             :   type, public :: hist_coord_t
         268             :     character(len=max_hcoordname_len) :: name = ''  ! coordinate name
         269             :     integer                  :: dimsize = 0       ! size of dimension
         270             :     character(len=max_hcoordname_len) :: dimname = '' ! optional dimension name
         271             :     character(len=max_chars) :: long_name = ''    ! 'long_name' attribute
         272             :     character(len=max_chars) :: units = ''        ! 'units' attribute
         273             :     character(len=max_chars) :: bounds_name = ''  ! 'bounds' attribute (& name of bounds variable)
         274             :     character(len=max_chars) :: standard_name = ''! 'standard_name' attribute
         275             :     character(len=4)         :: positive = ''     ! 'positive' attribute ('up' or 'down')
         276             :     integer,  pointer        :: integer_values(:) => null() ! dim values if integral
         277             :     real(r8), pointer        :: real_values(:) => null() ! dim values if real
         278             :     real(r8), pointer        :: bounds(:,:) => null() ! dim bounds
         279             :     type(formula_terms_t)    :: formula_terms     ! vars for formula terms
         280             :     logical                  :: integer_dim       ! .true. iff dim has integral values
         281             :     logical                  :: vertical_coord    ! .true. iff dim is vertical
         282             :   end type hist_coord_t
         283             :
         284             :   ! Some parameters for use with interpolated output namelist items
         285             :   integer,          parameter, public :: interp_type_native            = 0
         286             :   integer,          parameter, public :: interp_type_bilinear          = 1
         287             :   integer,          parameter, public :: interp_gridtype_equal_poles   = 1
         288             :   integer,          parameter, public :: interp_gridtype_gauss         = 2
         289             :   integer,          parameter, public :: interp_gridtype_equal_nopoles = 3
         290             :
         291             :   !---------------------------------------------------------------------------
         292             :   !
         293             :   !  interp_info_t: Information for lat/lon interpolated history output
         294             :   !
         295             :   !---------------------------------------------------------------------------
         296             :   type, public :: interp_info_t
         297             :     ! store the  lat-lon grid information
         298             :     character(len=28)     :: gridname = ''
         299             :     integer               :: grid_id  = -1
         300             :     ! gridtype = 1      equally spaced, including poles (FV scalars output grid)
         301             :     ! gridtype = 2      Gauss grid (CAM Eulerian)
         302             :     ! gridtype = 3      equally spaced, no poles (FV staggered velocity)
         303             :     integer               :: interp_gridtype = interp_gridtype_equal_poles
         304             :     ! interpolate_type = 0: native high order interpolation
         305             :     ! interpolate_type = 1: bilinear interpolation
         306             :     integer               :: interp_type = interp_type_bilinear
         307             :     integer               :: interp_nlat = 0
         308             :     integer               :: interp_nlon = 0
         309             :     real(r8), pointer     :: interp_lat(:) => NULL()
         310             :     real(r8), pointer     :: interp_lon(:) => NULL()
         311             :     real(r8), pointer     :: interp_gweight(:) => NULL()
         312             :   end type interp_info_t
         313             :
         314             :   !! Coordinate variables
         315             :   integer,                     public :: registeredmdims = 0
         316             :   integer,                     public :: maxvarmdims     = 1
         317             :   character(len=9), parameter, public :: mdim_var_name   = 'mdimnames'
         318             :   integer,          parameter         :: maxmdims        = 25  ! arbitrary limit
         319             :   type(hist_coord_t),          public :: hist_coords(maxmdims)
         320             :
         321             :   public     :: add_hist_coord, add_vert_coord
         322             :   public     :: write_hist_coord_attrs, write_hist_coord_vars
         323             :   public     :: lookup_hist_coord_indices, hist_coord_find_levels
         324             :   public     :: get_hist_coord_index, hist_coord_name, hist_coord_size
         325             :   public     :: hist_dimension_name
         326             :   public     :: hist_dimension_values
         327             :
         328             :   interface add_hist_coord
         329             :     module procedure add_hist_coord_regonly
         330             :     module procedure add_hist_coord_int
         331             :     module procedure add_hist_coord_r8
         332             :   end interface
         333             :
         334             :   interface hist_coord_size
         335             :     module procedure hist_coord_size_char
         336             :     module procedure hist_coord_size_int
         337             :  end interface hist_coord_size
         338             :
         339             :  interface hist_dimension_values
         340             :     module procedure hist_dimension_values_r8
         341             :     module procedure hist_dimension_values_int
         342             :  end interface hist_dimension_values
         343             :
         344             :   interface assignment(=)
         345             :     module procedure field_copy
         346             :     module procedure formula_terms_copy
         347             :   end interface
         348             :
         349             :   interface check_hist_coord
         350             :     ! NB: This is supposed to be a private interface
         351             :     ! check_hist_coord: returns 0 if <name> is not registered as an mdim
         352             :     !                   returns i if <name> is registered with compatible values
         353             :     !              calls endrun if <name> is registered with incompatible values
         354             :     ! Versions without the <name> argument return .true. or .false.
         355             :     module procedure check_hist_coord_char
         356             :     module procedure check_hist_coord_int
         357             :     module procedure check_hist_coord_int1
         358             :     module procedure check_hist_coord_r8
         359             :     module procedure check_hist_coord_r81
         360             :     module procedure check_hist_coord_r82
         361             :     module procedure check_hist_coord_ft
         362             :     module procedure check_hist_coord_all
         363             :   end interface
         364             :
         365             : !!---------------------------------------------------------------------------
         366             :
         367             : contains
         368             :
         369   102353760 :   subroutine dim_index_2d_dim_sizes_2d(this, dim1, dim2)
         370             :
         371             :     ! Dummy arguments
         372             :     class(dim_index_2d)                :: this
         373             :     integer,             intent(out)   :: dim1
         374             :     integer,             intent(out)   :: dim2
         375             :
         376   102353760 :     dim1 = MAX(0, this%end1 - this%beg1 + 1)
         377   102353760 :     dim2 = MAX(0, this%end2 - this%beg2 + 1)
         378             :
         379   102353760 :   end subroutine dim_index_2d_dim_sizes_2d
         380             :
         381           0 :   subroutine dim_index_2d_dim_size_arr(this, dims)
         382             :
         383             :     ! Dummy arguments
         384             :     class(dim_index_2d)                :: this
         385             :     integer,             intent(out)   :: dims(:)
         386             :
         387           0 :     if (size(dims) < 2) then
         388           0 :       call endrun('dim_index_2d_dim_sizes: dims must have at least two elements')
         389             :     end if
         390             :
         391           0 :     call this%dim_sizes(dims(1), dims(2))
         392             :
         393           0 :   end subroutine dim_index_2d_dim_size_arr
         394             :
         395    10690560 :   subroutine dim_index_3d_dim_sizes_3d(this, dim1, dim2, dim3)
         396             :
         397             :     ! Dummy arguments
         398             :     class(dim_index_3d)                :: this
         399             :     integer,             intent(out)   :: dim1
         400             :     integer,             intent(out)   :: dim2
         401             :     integer,             intent(out)   :: dim3
         402             :
         403    10690560 :     dim1 = MAX(0, this%end1 - this%beg1 + 1)
         404    10690560 :     dim2 = MAX(0, this%end2 - this%beg2 + 1)
         405    10690560 :     dim3 = MAX(0, this%end3 - this%beg3 + 1)
         406             :
         407    10690560 :   end subroutine dim_index_3d_dim_sizes_3d
         408             :
         409    10690560 :   subroutine dim_index_3d_dim_size_arr(this, dims)
         410             :
         411             :     ! Dummy arguments
         412             :     class(dim_index_3d)                :: this
         413             :     integer,             intent(out)   :: dims(:)
         414             :
         415    10690560 :     if (size(dims) < 3) then
         416           0 :       call endrun('dim_index_3d_dim_sizes: dims must have at least three elements')
         417             :     end if
         418             :
         419    10690560 :     call this%dim_sizes(dims(1), dims(2), dims(3))
         420             :
         421    10690560 :   end subroutine dim_index_3d_dim_size_arr
         422             :
         423             :   ! field_info_get_dims_2d: Retrieve bounds for stepping through a chunk
         424   231723864 :   type(dim_index_2d) function field_info_get_dims_2d(this, col) result(dims)
         425             :     use cam_grid_support, only: cam_grid_get_block_count
         426             :
         427             :     ! Dummy argument
         428             :     class(field_info)                :: this
         429             :     integer,           intent(in)    :: col
         430             :
         431             :     ! Local variable
         432             :     integer              :: endi
         433             :
         434   231723864 :     if (this%colperchunk) then
         435   231686064 :       endi = this%begdim1 + cam_grid_get_block_count(this%decomp_type, col) - 1
         436   231686064 :       dims = dim_index_2d(this%begdim1, endi, this%begdim2, this%enddim2)
         437             :     else
         438       37800 :       dims = dim_index_2d(this%begdim1, this%enddim1, this%begdim2, this%enddim2)
         439             :     end if
         440   231723864 :   end function field_info_get_dims_2d
         441             :
         442             :   ! field_info_get_dims_3d: Retrieve grid decomp bounds
         443    10690560 :   type(dim_index_3d) function field_info_get_dims_3d(this) result(dims)
         444             :
         445             :     ! Dummy argument
         446             :     class(field_info)                :: this
         447             :
         448             :     dims = dim_index_3d(this%begdim1, this%enddim1, this%begdim2, this%enddim2,&
         449    10690560 :          this%begdim3, this%enddim3)
         450             :
         451   231723864 :   end function field_info_get_dims_3d
         452             :
         453             :   ! field_info_is_composed: Return whether this field is composed of two other fields
         454    21381120 :   pure logical function field_info_is_composed(this)
         455             :     class(field_info), intent(IN)  :: this
         456             :
         457             :     field_info_is_composed = ((trim(adjustl(this%field_op))=='sum' .or. trim(adjustl(this%field_op))=='dif') .and. &
         458    21381120 :                               this%op_field1_id /= -1 .and. this%op_field2_id /= -1)
         459    21381120 :   end function field_info_is_composed
         460             :
         461             :   ! field_info_get_shape: Return a pointer to the field's global shape.
         462             :   !                       Calculate it first if necessary
         463    10690560 :   subroutine field_info_get_shape(this, shape_out, rank_out)
         464             :     use cam_grid_support, only: cam_grid_dimensions
         465             :
         466             :     ! Dummy arguments
         467             :     class(field_info)                         :: this
         468             :     integer,                    intent(out)   :: shape_out(:)
         469             :     integer,                    intent(out)   :: rank_out
         470             :
         471             :     ! Local arguments
         472             :     integer                                   :: rank, i, pos
         473             :     integer                                   :: gdims(2)
         474             :
         475    10690560 :     if (.not. associated(this%shape)) then
         476             :       ! Calculate field's global shape
         477      133632 :       call cam_grid_dimensions(this%decomp_type, gdims, rank)
         478      133632 :       pos = rank
         479      133632 :       if (associated(this%mdims)) then
         480       88320 :         rank = rank + size(this%mdims)
         481             :       end if
         482      400896 :       allocate(this%shape(rank))
         483      267264 :       this%shape(1:pos) = gdims(1:pos)
         484      267264 :       if (rank > pos) then
         485       86016 :         do i = 1, size(this%mdims)
         486       43008 :           pos = pos + 1
         487       86016 :           this%shape(pos) = hist_coords(this%mdims(i))%dimsize
         488             :         end do
         489             :       end if
         490             :     end if
         491             :
         492    10690560 :     rank_out = size(this%shape)
         493             :
         494    10690560 :     if (size(shape_out) < rank_out) then
         495           0 :       call endrun('field_info_get_shape: shape_out too small')
         496             :     end if
         497             :
         498    24821760 :     shape_out(1:rank_out) = this%shape(1:rank_out)
         499    10690560 :     if (size(shape_out) > rank_out) then
         500    92774400 :       shape_out(rank_out+1:) = 1
         501             :     end if
         502             :
         503    10690560 :   end subroutine field_info_get_shape
         504             :
         505    32087808 :   subroutine field_info_get_bounds(this, dim, beg, end)
         506             :
         507             :     ! Dummy arguments
         508             :     class(field_info)                :: this
         509             :     integer,           intent(in)    :: dim
         510             :     integer,           intent(out)   :: beg
         511             :     integer,           intent(out)   :: end
         512             :
         513    32087808 :     select case(dim)
         514             :     case (1)
         515           0 :       beg = this%begdim1
         516           0 :       end = this%enddim1
         517             :     case (2)
         518           0 :       beg = this%begdim2
         519           0 :       end = this%enddim2
         520             :     case (3)
         521    32087808 :       beg = this%begdim3
         522    32087808 :       end = this%enddim3
         523             :     case default
         524    32087808 :       call endrun('field_info_get_bounds: dim must be 1, 2, or 3')
         525             :     end select
         526             :
         527    10690560 :   end subroutine field_info_get_bounds
         528             :
         529      144384 :   subroutine hentry_get_global(this, gval)
         530             :
         531             :     ! Dummy arguments
         532             :     class(hentry)                    :: this
         533             :     real(r8),          intent(out)   :: gval
         534             :
         535      144384 :     gval=this%hbuf_integral
         536             :
         537      144384 :   end subroutine hentry_get_global
         538             :
         539       72192 :   subroutine hentry_put_global(this, gval)
         540             :
         541             :     ! Dummy arguments
         542             :     class(hentry)                    :: this
         543             :     real(r8),          intent(in)    :: gval
         544             :
         545       72192 :     this%hbuf_integral=gval
         546             :
         547       72192 :   end subroutine hentry_put_global
         548             :
         549             :   ! history_patch_write_attrs: Define coordinate variables and attributes
         550             :   !               for a patch
         551           0 :   subroutine history_patch_write_attrs(this, File)
         552             :     use cam_grid_support, only: cam_grid_is_unstructured
         553             :     use pio,              only: file_desc_t, var_desc_t, pio_put_att, pio_double
         554             :     use cam_pio_utils,    only: cam_pio_def_dim, cam_pio_def_var, cam_pio_handle_error
         555             :
         556             :     ! Dummy arguments
         557             :     class(history_patch_t)                  :: this
         558             :     type(file_desc_t),        intent(inout) :: File    ! PIO file Handle
         559             :
         560             :     ! Local variable
         561             :     type(cam_grid_patch_t), pointer         :: patchptr
         562             :     type(var_desc_t), pointer               :: vardesc_lat => NULL()
         563             :     type(var_desc_t), pointer               :: vardesc_lon => NULL()
         564             :     character(len=128)                      :: errormsg
         565             :     character(len=max_chars)                :: lat_name
         566             :     character(len=max_chars)                :: lon_name
         567             :     character(len=max_chars)                :: col_name
         568             :     character(len=max_chars)                :: temp_str
         569             :     integer                                 :: dimid1, dimid2 ! PIO dim IDs
         570             :     integer                                 :: num_patches
         571             :     integer                                 :: temp1, temp2
         572             :     integer                                 :: latid, lonid ! Coordinate dims
         573             :     integer                                 :: i, ierr
         574             :     logical                                 :: col_only
         575             :     logical                                 :: unstruct
         576             :     character(len=*), parameter             :: subname = 'history_patch_write_attrs'
         577             :
         578           0 :     num_patches = size(this%patches)
         579           0 :     if (associated(this%header_info)) then
         580             :       ! Make sure header_info is the right size
         581           0 :       if (size(this%header_info) /= num_patches) then
         582           0 :         write(errormsg, '(a,2(i0,a))') 'Size mismatch between header_info (', &
         583           0 :              size(this%header_info), ') and patches (', num_patches, ')'
         584           0 :         call endrun(subname//': '//errormsg)
         585             :       end if
         586             :     else
         587           0 :       allocate(this%header_info(num_patches))
         588             :     end if
         589             :
         590             :     ! Write attributes for each patch
         591           0 :     do i = 1, num_patches
         592           0 :       patchptr => this%patches(i)
         593           0 :       call this%header_info(i)%set_gridid(patchptr%gridid())
         594           0 :       unstruct = cam_grid_is_unstructured(patchptr%gridid())
         595             :       ! What are the dimension(s) for this patch?
         596           0 :       col_only = this%collected_output
         597           0 :       if (num_patches == 1) then
         598             :         ! Backwards compatibility
         599           0 :         if (unstruct .or. col_only) then
         600           0 :           col_name = 'ncol'
         601             :         else
         602           0 :           col_name = ''
         603             :         end if
         604           0 :         lat_name = 'lat'
         605           0 :         lon_name = 'lon'
         606             :       else
         607           0 :         call patchptr%get_axis_names(lat_name, lon_name, col_name, col_only)
         608             :       end if
         609             :       ! Define the dimensions (latx/lonx or ncolx)
         610             :       ! col_name is set for unstructured output or collected columns (ncolx)
         611           0 :       if (len_trim(col_name) > 0) then
         612           0 :         call patchptr%get_global_size(gsize=temp1)
         613           0 :         if (temp1 <= 0) then
         614           0 :           call endrun(subname//': col dimsize must be positive')
         615             :         end if
         616           0 :         if (unstruct .and. (.not. col_only)) then
         617             :           ! For the case of unstructured output without collected column
         618             :           ! output, we need to make sure that the ncolx dimension is unique
         619           0 :           col_name = trim(col_name)//'_'//trim(this%lon_axis_name)//'_'//trim(this%lat_axis_name)
         620             :         end if
         621           0 :         call cam_pio_def_dim(File, trim(col_name), temp1, dimid1, existOK=.false.)
         622           0 :         call this%header_info(i)%set_hdims(dimid1)
         623           0 :         latid = dimid1
         624           0 :         lonid = dimid1
         625             :       else
         626           0 :         lat_name = trim(lat_name)//'_'//trim(this%lat_axis_name)
         627           0 :         call patchptr%get_global_size(temp1, temp2)
         628           0 :         if (temp1 <= 0) then
         629           0 :           call endrun(subname//': lat dimsize must be positive')
         630             :         end if
         631           0 :         call cam_pio_def_dim(File, trim(lat_name), temp1, dimid1, existOK=.true.)
         632           0 :         latid = dimid1
         633           0 :         lon_name = trim(lon_name)//'_'//trim(this%lon_axis_name)
         634           0 :         if (temp2 <= 0) then
         635           0 :           call endrun(subname//': lon dimsize must be positive')
         636             :         end if
         637           0 :         call cam_pio_def_dim(File, trim(lon_name), temp2, dimid2, existOK=.true.)
         638           0 :         lonid = dimid2
         639           0 :         call this%header_info(i)%set_hdims(lonid, latid)
         640             :       end if
         641             :       !! Define the latx (coordinate) variable
         642           0 :       if (unstruct .and. (.not. col_only)) then
         643             :         ! We need to make sure the latx name is unique
         644           0 :         lat_name = trim(lat_name)//'_'//trim(this%lon_axis_name)//'_'//trim(this%lat_axis_name)
         645             :       end if
         646           0 :       allocate(vardesc_lat)
         647             :       call cam_pio_def_var(File, trim(lat_name), pio_double, (/latid/),       &
         648           0 :            vardesc_lat, existOK=.true.)
         649             :       ! Coordinate attributes
         650           0 :       call patchptr%get_coord_long_name('lat', temp_str)
         651           0 :       if (len_trim(temp_str) > 0) then
         652           0 :         ierr = pio_put_att(File, vardesc_lat, 'long_name', trim(temp_str))
         653           0 :         call cam_pio_handle_error(ierr, subname//': Unable to define long_name')
         654             :       end if
         655           0 :       call patchptr%get_coord_units('lat', temp_str)
         656           0 :       if (len_trim(temp_str) > 0) then
         657           0 :         ierr = pio_put_att(File, vardesc_lat, 'units', trim(temp_str))
         658           0 :         call cam_pio_handle_error(ierr, subname//': Unable to define units')
         659             :       end if
         660             :       !! Define the lonx (coordinate) variable
         661           0 :       if (unstruct .and. (.not. col_only)) then
         662             :         ! We need to make sure the lonx name is unique
         663           0 :         lon_name = trim(lon_name)//'_'//trim(this%lon_axis_name)//'_'//trim(this%lat_axis_name)
         664             :       end if
         665           0 :       allocate(vardesc_lon)
         666             :       call cam_pio_def_var(File, trim(lon_name), pio_double, (/lonid/),       &
         667           0 :            vardesc_lon, existOK=.true.)
         668             :       ! Coordinate attributes
         669           0 :       call patchptr%get_coord_long_name('lon', temp_str)
         670           0 :       if (len_trim(temp_str) > 0) then
         671           0 :         ierr = pio_put_att(File, vardesc_lon, 'long_name', trim(temp_str))
         672           0 :         call cam_pio_handle_error(ierr, subname//': Unable to define long_name')
         673             :       end if
         674           0 :       call patchptr%get_coord_units('lon', temp_str)
         675           0 :       if (len_trim(temp_str) > 0) then
         676           0 :         ierr = pio_put_att(File, vardesc_lon, 'units', trim(temp_str))
         677           0 :         call cam_pio_handle_error(ierr, subname//': Unable to define units')
         678             :       end if
         679           0 :       call this%header_info(i)%set_varids(vardesc_lon, vardesc_lat)
         680           0 :       nullify(vardesc_lat, vardesc_lon) ! They belong to the header_info now
         681             :     end do
         682             :
         683           0 :   end subroutine history_patch_write_attrs
         684             :
         685             :   ! history_patch_write_vals: Write coordinate variable values for a patch
         686           0 :   subroutine history_patch_write_vals(this, File)
         687           0 :     use pio,           only: file_desc_t, var_desc_t
         688             :
         689             :     ! Dummy arguments
         690             :     class(history_patch_t)                  :: this
         691             :     type(file_desc_t),        intent(inout) :: File    ! PIO file Handle
         692             :
         693             :     ! Local variable
         694             :     type(cam_grid_patch_t), pointer         :: patchptr
         695             :     type(var_desc_t), pointer               :: vardesc => NULL()  ! PIO var desc
         696             :     character(len=128)                      :: errormsg
         697             :     integer                                 :: num_patches
         698             :     integer                                 :: i
         699             :
         700           0 :     num_patches = size(this%patches)
         701           0 :     if (.not. associated(this%header_info)) then
         702             :       ! We need this for dim and coord variable descriptors
         703           0 :       write(errormsg, '(2a)') 'No header info for ', trim(this%namelist_entry)
         704           0 :       call endrun('history_patch_write_vals: '//errormsg)
         705             :     end if
         706             :
         707             :     ! Write attributes for each patch
         708           0 :     do i = 1, num_patches
         709           0 :       patchptr => this%patches(i)
         710             :       ! Write the coordinate variables (or just lat/lon for column output)
         711           0 :       call patchptr%write_coord_vals(File, this%header_info(i))
         712             :     end do
         713             :
         714           0 :   end subroutine history_patch_write_vals
         715             :
         716             :   ! history_patch_field_name: Add patch description to field name
         717           0 :   subroutine history_patch_field_name(this, name)
         718             :     ! Dummy arguments
         719             :     class(history_patch_t)                     :: this
         720             :     character(len=*),            intent(inout) :: name
         721             :
         722           0 :     if (.not. this%collected_output) then
         723             :       ! Add patch description info to the variable name
         724           0 :       name = trim(name)//'_'//trim(this%lon_axis_name)//'_'//trim(this%lat_axis_name)
         725             :     end if
         726           0 :   end subroutine history_patch_field_name
         727             :
         728             :   ! history_patch_num_hdims: Find the number of horizontal dimensions for
         729             :   !         the indicated grid
         730           0 :   integer function history_patch_num_hdims(this, gridid)
         731             :     ! Dummy arguments
         732             :     class(history_patch_t)                     :: this
         733             :     integer,                     intent(in)    :: gridid      ! The field's grid
         734             :
         735             :     ! Local variables
         736             :     type(cam_grid_patch_t),        pointer     :: patchptr
         737             :     character(len=128)                         :: errormsg
         738             :     integer                                    :: i
         739             :     integer                                    :: num_patches
         740             :
         741             :     ! Basic sanity checks, is this patch OK?
         742           0 :     num_patches = size(this%patches)
         743           0 :     if (associated(this%header_info)) then
         744             :       ! Make sure header_info is the right size
         745           0 :       if (size(this%header_info) /= num_patches) then
         746           0 :         write(errormsg, '(a,2(i0,a))') 'Size mismatch between header_info (', &
         747           0 :              size(this%header_info), ') and patches (', num_patches, ')'
         748           0 :         call endrun('history_patch_num_hdims: '//errormsg)
         749             :       end if
         750             :     else
         751           0 :       write(errormsg, *) 'No header info for patch, ', trim(this%namelist_entry)
         752           0 :       call endrun('history_patch_num_hdims: '//errormsg)
         753             :     end if
         754             :
         755             :     ! Find the correct patch by matching grid ID
         756           0 :     history_patch_num_hdims = -1
         757           0 :     do i = 1, num_patches
         758           0 :       patchptr => this%patches(i)
         759           0 :       if (patchptr%gridid() == gridid) then
         760             :         ! This is the right patch, set the return val and quit loop
         761           0 :         history_patch_num_hdims = this%header_info(i)%num_hdims()
         762           0 :         exit
         763           0 :       else if (i >= num_patches) then
         764           0 :         write(errormsg, '(3a,i0)') 'No grid found for patch, ',               &
         765           0 :              trim(this%namelist_entry), '. Was looking for decomp ', gridid
         766           0 :         call endrun('history_patch_num_hdims: '//errormsg)
         767             :       ! No else needed
         768             :       end if
         769             :     end do
         770           0 :     if (history_patch_num_hdims <= 0) then
         771           0 :       write(errormsg, '(2a,2(a,i0))') 'INTERNAL: No grid patch for ',         &
         772           0 :            trim(this%namelist_entry), ', num_patches = ',num_patches,         &
         773           0 :            ', gridid = ', gridid
         774           0 :       call endrun('history_patch_num_hdims: '//errormsg)
         775             :     end if
         776             :
         777           0 :   end function history_patch_num_hdims
         778             :
         779             :   ! history_patch_get_var_data: Calculate data relevant to history variable
         780             :   !         on a patch by substituting patch dimension ids for the horiz. ids
         781             :   !         and adding patch information to the variable name
         782           0 :   subroutine history_patch_get_var_data(this, name, dimids, gridid)
         783             :     ! Dummy arguments
         784             :     class(history_patch_t)                     :: this
         785             :     character(len=*),            intent(inout) :: name
         786             :     integer,                     intent(inout) :: dimids(:)   ! Grid dimids
         787             :     integer,                     intent(in)    :: gridid      ! The field's grid
         788             :
         789             :     ! Local variables
         790             :     type(cam_grid_patch_t),        pointer     :: patchptr
         791             :     type (cam_grid_header_info_t), pointer     :: histptr
         792             :     character(len=128)                         :: errormsg
         793             :     integer                                    :: num_patches
         794             :     integer                                    :: i
         795             :
         796             :     ! Basic sanity checks, is this patch OK?
         797           0 :     num_patches = size(this%patches)
         798           0 :     if (associated(this%header_info)) then
         799             :       ! Make sure header_info is the right size
         800           0 :       if (size(this%header_info) /= num_patches) then
         801           0 :         write(errormsg, '(a,2(i0,a))') 'Size mismatch between header_info (', &
         802           0 :              size(this%header_info), ') and patches (', num_patches, ')'
         803           0 :         call endrun('history_patch_get_var_data: '//errormsg)
         804             :       end if
         805             :     else
         806           0 :       write(errormsg, *) 'No header info for patch, ', trim(this%namelist_entry)
         807           0 :       call endrun('history_patch_get_var_data: '//errormsg)
         808             :     end if
         809             :
         810             :     ! Find the correct patch by matching grid ID
         811           0 :     do i = 1, num_patches
         812           0 :       patchptr => this%patches(i)
         813           0 :       if (patchptr%gridid() == gridid) then
         814             :         ! This is the right patch, quit loop
         815           0 :         histptr  => this%header_info(i)
         816           0 :         exit
         817           0 :       else if (i >= num_patches) then
         818           0 :         write(errormsg, '(3a,i0)') 'No grid found for patch, ',               &
         819           0 :              trim(this%namelist_entry), '. Was looking for decomp ', gridid
         820           0 :         call endrun('history_patch_get_var_data: '//errormsg)
         821             :       ! No else needed
         822             :       end if
         823             :     end do
         824             :
         825             :     ! We have the correct patch, replace the horizontal dimension ids
         826           0 :     do i = 1, histptr%num_hdims()
         827           0 :       dimids(i) = histptr%get_hdimid(i)
         828             :     end do
         829             :     ! Re-define the variable name
         830           0 :     call this%field_name(name)
         831             :
         832           0 :   end subroutine history_patch_get_var_data
         833             :
         834           0 :   subroutine history_patch_compact(this)
         835             :
         836             :     ! Dummy arguments
         837             :     class(history_patch_t)                   :: this
         838             :
         839             :     ! Local variables
         840             :     integer                                  :: num_patches
         841             :     integer                                  :: i
         842             :
         843           0 :     num_patches = size(this%patches)
         844             :
         845             :     ! Find the correct patch by matching grid ID
         846           0 :     do i = 1, num_patches
         847           0 :       call this%patches(i)%compact(this%collected_output)
         848             :     end do
         849             :
         850           0 :   end subroutine history_patch_compact
         851             :
         852           0 :   subroutine history_patch_write_var(this, File, gridid, adims, dtype, hbuf, varid)
         853             :     use pio,           only: file_desc_t, var_desc_t, io_desc_t
         854             :     use pio,           only: pio_write_darray
         855             :     use cam_pio_utils, only: cam_pio_handle_error, cam_pio_var_info
         856             :
         857             :     ! Dummy arguments
         858             :     class(history_patch_t)                   :: this
         859             :     type(file_desc_t),         intent(inout) :: File       ! PIO file handle
         860             :     integer,                   intent(in)    :: gridid
         861             :     integer,                   intent(in)    :: adims(:)
         862             :     integer,                   intent(in)    :: dtype
         863             :     real(r8),                  intent(in)    :: hbuf(:,:,:)
         864             :     type(var_desc_t),       pointer          :: varid
         865             :
         866             :     ! Local variables
         867             :     type(cam_grid_patch_t), pointer          :: patchptr
         868             :     character(len=128)                       :: errormsg
         869             :     integer                                  :: num_patches
         870             :     integer                                  :: i, idx
         871             :     integer                                  :: uid ! unlimited dim ID
         872             :     type(io_desc_t),        pointer          :: iodesc
         873             :     integer                                  :: ierr, nfdims
         874             :     integer                                  :: fdimlens(7), dimids(7)
         875             :
         876           0 :     num_patches = size(this%patches)
         877             :
         878             :     ! Find the correct patch by matching grid ID
         879           0 :     do i = 1, num_patches
         880           0 :       patchptr => this%patches(i)
         881           0 :       if (patchptr%gridid() == gridid) then
         882             :         ! This is the right patch, quit loop
         883             :         exit
         884           0 :       else if (i >= num_patches) then
         885           0 :         write(errormsg, '(3a,i0)') 'No grid found for patch, ',               &
         886           0 :              trim(this%namelist_entry), '. Was looking for decomp ', gridid
         887           0 :         call endrun('history_patch_write_var: '//trim(errormsg))
         888             :       ! No else needed
         889             :       end if
         890             :     end do
         891             :
         892             :     ! We have the right grid, write the hbuf
         893           0 :     call cam_pio_var_info(File, varid, nfdims, dimids, fdimlens, unlimDimID=uid)
         894           0 :     idx = 1
         895           0 :     do i = 1, nfdims
         896           0 :       if (i > idx) then
         897           0 :         dimids(idx) = dimids(i)
         898             :       end if
         899           0 :       if (dimids(i) /= uid) then
         900           0 :         idx = idx + 1
         901             :       end if
         902             :     end do
         903           0 :     nfdims = nfdims - COUNT(dimids(1:nfdims) == uid)
         904           0 :     call patchptr%get_decomp(adims, fdimlens(1:nfdims), dtype, iodesc)
         905           0 :     if (size(adims) == 2) then
         906           0 :       call pio_write_darray(File, varid, iodesc, hbuf(:,1,:), ierr)
         907           0 :     else if (size(adims) == 3) then
         908           0 :       call pio_write_darray(File, varid, iodesc, hbuf, ierr)
         909             :     else
         910           0 :       call endrun("history_patch_write_var: adims must be rank 2 or 3")
         911             :     end if
         912           0 :     call cam_pio_handle_error(ierr, 'history_patch_write_var: Error writing variable')
         913             :
         914           0 :   end subroutine history_patch_write_var
         915             :
         916           0 :   subroutine history_patch_active_cols(this, gridid, lchnk, active)
         917             :     ! Dummy arguments
         918             :     class(history_patch_t)                :: this
         919             :     integer,                intent(in)    :: gridid ! desired grid
         920             :     integer,                intent(in)    :: lchnk  ! chunk or block number
         921             :     logical,                intent(out)   :: active(:)
         922             :
         923             :     ! Local variables
         924             :     type(cam_grid_patch_t),        pointer     :: patchptr
         925             :     character(len=128)                         :: errormsg
         926             :     integer                                    :: num_patches
         927             :     integer                                    :: i
         928             :
         929           0 :     num_patches = size(this%patches)
         930             :
         931             :     ! Find the correct patch by matching grid ID
         932           0 :     do i = 1, num_patches
         933           0 :       patchptr => this%patches(i)
         934           0 :       if (patchptr%gridid() == gridid) then
         935             :         ! This is the right patch, quit loop
         936             :         exit
         937           0 :       else if (i >= num_patches) then
         938           0 :         write(errormsg, '(3a,i0)') 'No grid found for patch, ',               &
         939           0 :              trim(this%namelist_entry), '. Was looking for decomp ', gridid
         940           0 :         call endrun('history_patch_active_cols: '//errormsg)
         941             :       ! No else needed
         942             :       end if
         943             :     end do
         944             :
         945             :     ! If we get here, patchptr is the grid patch we want
         946           0 :     call patchptr%active_cols(lchnk, active)
         947             :
         948           0 :  end subroutine history_patch_active_cols
         949             :
         950           0 :   subroutine history_patch_deallocate(this)
         951             :     ! Dummy argument
         952             :     class(history_patch_t) :: this
         953             :     ! Local variable
         954             :     integer                :: i
         955             :
         956           0 :     this%lon_axis_name = ''
         957           0 :     this%lat_axis_name = ''
         958             :
         959           0 :     if (associated(this%patches)) then
         960           0 :       do i = 1, size(this%patches)
         961           0 :         call this%patches(i)%deallocate()
         962             :       end do
         963           0 :       deallocate(this%patches)
         964           0 :       nullify(this%patches)
         965             :     end if
         966             :
         967           0 :     if (associated(this%header_info)) then
         968           0 :       do i = 1, size(this%header_info)
         969           0 :         call this%header_info(i)%deallocate()
         970             :       end do
         971           0 :       deallocate(this%header_info)
         972           0 :       nullify(this%header_info)
         973             :     end if
         974             :
         975           0 :   end subroutine history_patch_deallocate
         976             :
         977             :   subroutine field_copy(f_out, f_in)
         978             :     type(field_info), intent(in) :: f_in
         979             :     type(field_info), intent(out) :: f_out
         980             :
         981             :     f_out%flag_xyfill= f_in%flag_xyfill
         982             :     f_out%is_subcol = f_in%is_subcol
         983             :     f_out%fillvalue= f_in%fillvalue
         984             :     f_out%numlev =  f_in%numlev                      ! vertical dimension (.nc file and internal arr)
         985             :     f_out%begdim1 = f_in%begdim1                     ! on-node dim1 start index
         986             :     f_out%enddim1 = f_in%enddim1                     ! on-node dim1 end index
         987             :     f_out%begdim2 = f_in%begdim2                     ! on-node dim2 start index
         988             :     f_out%enddim2 = f_in%enddim2                     ! on-node dim2 end index
         989             :     f_out%begdim3 = f_in%begdim3                     ! on-node chunk or lat start index
         990             :     f_out%enddim3 = f_in%enddim3                     ! on-node chunk or lat end index
         991             :     f_out%decomp_type = f_in%decomp_type             ! type of decomposition (physics or dynamics)
         992             :
         993             :     f_out%meridional_complement = f_in%meridional_complement ! id  or -1
         994             :     f_out%zonal_complement = f_in%zonal_complement           ! id  or -1
         995             :     f_out%field_op = f_in%field_op                           ! sum,dif, or ''
         996             :     f_out%op_field1_id = f_in%op_field1_id   ! id  or -1
         997             :     f_out%op_field2_id = f_in%op_field2_id   ! id  or -1
         998             :
         999             :     f_out%name = f_in%name                           ! field name
        1000             :     f_out%long_name = f_in%long_name                 ! long name
        1001             :     f_out%units = f_in%units                         ! units
        1002             :     f_out%mixing_ratio = f_in%mixing_ratio           ! mixing_ratio
        1003             :     f_out%sampling_seq =  f_in%sampling_seq          ! sampling sequence - if not every timestep, how often field is sampled
        1004             :     f_out%cell_methods = f_in%cell_methods
        1005             :
        1006             :     if(associated(f_in%mdims)) then
        1007             :       f_out%mdims=>f_in%mdims
        1008             :     else
        1009             :       nullify(f_out%mdims)
        1010             :     end if
        1011             :
        1012             :   end subroutine field_copy
        1013             :
        1014           0 :   subroutine formula_terms_copy(f_out, f_in)
        1015             :     type(formula_terms_t), intent(in) :: f_in
        1016             :     type(formula_terms_t), intent(out) :: f_out
        1017             :
        1018           0 :     f_out%a_name = f_in%a_name
        1019           0 :     f_out%a_long_name = f_in%a_long_name
        1020           0 :     f_out%a_values => f_in%a_values
        1021           0 :     f_out%b_name = f_in%b_name
        1022           0 :     f_out%b_long_name = f_in%b_long_name
        1023           0 :     f_out%b_values => f_in%b_values
        1024           0 :     f_out%p0_name = f_in%p0_name
        1025           0 :     f_out%p0_long_name = f_in%p0_long_name
        1026           0 :     f_out%p0_units = f_in%p0_units
        1027           0 :     f_out%p0_value = f_in%p0_value
        1028           0 :     f_out%ps_name = f_in%ps_name
        1029           0 :   end subroutine formula_terms_copy
        1030             :
        1031      683520 :   integer function get_hist_coord_index(mdimname)
        1032             :     ! Input variables
        1033             :     character(len=*), intent(in)            :: mdimname
        1034             :     ! Local variable
        1035             :     integer :: i
        1036             :
        1037      683520 :     get_hist_coord_index = -1
        1038     1029120 :     do i = 1, registeredmdims
        1039     1029120 :       if(trim(mdimname) == trim(hist_coords(i)%name)) then
        1040      652800 :         get_hist_coord_index = i
        1041      652800 :         exit
        1042             :       end if
        1043             :     end do
        1044             :
        1045      683520 :   end function get_hist_coord_index
        1046             :
        1047       12288 :   character(len=max_hcoordname_len) function hist_coord_name(index)
        1048             :     ! Input variables
        1049             :     integer, intent(in)            :: index
        1050             :
        1051       12288 :     if ((index > 0) .and. (index <= registeredmdims)) then
        1052       12288 :       hist_coord_name = hist_coords(index)%name
        1053             :     else
        1054           0 :       call endrun('hist_coord_name: index out of range')
        1055             :     end if
        1056             :
        1057       12288 :   end function hist_coord_name
        1058             :
        1059      640512 :   integer function hist_coord_size_int(index)
        1060             :     ! Input variables
        1061             :     integer, intent(in)            :: index
        1062             :
        1063      640512 :     if (index > 0) then
        1064      640512 :       hist_coord_size_int = hist_coords(index)%dimsize
        1065             :     else
        1066             :       hist_coord_size_int = -1
        1067             :     end if
        1068             :
        1069      640512 :   end function hist_coord_size_int
        1070             :
        1071           0 :   integer function hist_coord_size_char(mdimname)
        1072             :     ! Input variables
        1073             :     character(len=*), intent(in)            :: mdimname
        1074             :     ! Local variable
        1075             :     integer :: i
        1076             :
        1077           0 :     i = get_hist_coord_index(mdimname)
        1078           0 :     hist_coord_size_char = hist_coord_size(i)
        1079             :
        1080           0 :   end function hist_coord_size_char
        1081             :
        1082             :   ! Functions to check consistent term definition for hist coords
        1083           0 :   logical function check_hist_coord_char(defined, input)
        1084             :
        1085             :     ! Input variables
        1086             :     character(len=*), intent(in)            :: defined
        1087             :     character(len=*), intent(in), optional  :: input
        1088             :
        1089           0 :     if (len_trim(defined) == 0) then
        1090             :       ! In this case, we assume the current value is undefined so any input OK
        1091             :       check_hist_coord_char = .true.
        1092           0 :     else if (present(input)) then
        1093             :       ! We have to match definitions
        1094           0 :       check_hist_coord_char = (trim(input) == trim(defined))
        1095             :     else
        1096             :       ! Not sure here. We have a value and are redefining without one?
        1097             :       check_hist_coord_char = .false.
        1098             :     end if
        1099           0 :   end function check_hist_coord_char
        1100             :
        1101           0 :   logical function check_hist_coord_int(defined, input)
        1102             :
        1103             :     ! Input variables
        1104             :     integer, intent(in)            :: defined
        1105             :     integer, intent(in), optional  :: input
        1106             :
        1107           0 :     if (defined == 0) then
        1108             :       ! In this case, we assume the current value is undefined so any input OK
        1109             :       check_hist_coord_int = .true.
        1110           0 :     else if (present(input)) then
        1111             :       ! We have to match definitions
        1112           0 :       check_hist_coord_int = (input == defined)
        1113             :     else
        1114             :       ! Not sure here. We have a value and are redefining without one?
        1115             :       check_hist_coord_int = .false.
        1116             :     end if
        1117           0 :   end function check_hist_coord_int
        1118             :
        1119           0 :   logical function check_hist_coord_int1(defined, input)
        1120             :
        1121             :     ! Input variables
        1122             :     integer,             pointer            :: defined(:)
        1123             :     integer, intent(in),          optional  :: input(:)
        1124             :
        1125             :     ! Local variables
        1126             :     integer                                 :: i
        1127             :
        1128           0 :     if (.not. associated(defined)) then
        1129             :       ! In this case, we assume the current value is undefined so any input OK
        1130             :       check_hist_coord_int1 = .true.
        1131           0 :     else if (present(input)) then
        1132             :       ! We have to match definitions
        1133           0 :       check_hist_coord_int1 = (size(input) == size(defined))
        1134             :     else
        1135             :       ! Not sure here. We have a value and are redefining without one?
        1136             :       check_hist_coord_int1 = .false.
        1137             :     end if
        1138           0 :     if (check_hist_coord_int1 .and. associated(defined)) then
        1139             :       ! Need to check the values
        1140           0 :       do i = 1, size(defined)
        1141           0 :         if (defined(i) /= input(i)) then
        1142             :           check_hist_coord_int1 = .false.
        1143             :           exit
        1144             :         end if
        1145             :       end do
        1146             :     end if
        1147           0 :   end function check_hist_coord_int1
        1148             :
        1149           0 :   logical function check_hist_coord_r8(defined, input)
        1150             :
        1151             :     ! Input variables
        1152             :     real(r8), intent(in)            :: defined
        1153             :     real(r8), intent(in), optional  :: input
        1154             :
        1155           0 :     if (defined == fillvalue) then
        1156             :       ! In this case, we assume the current value is undefined so any input OK
        1157             :       check_hist_coord_r8 = .true.
        1158           0 :     else if (present(input)) then
        1159             :       ! We have to match definitions
        1160           0 :       check_hist_coord_r8 = (input == defined)
        1161             :     else
        1162             :       ! Not sure here. We have a value and are redefining without one?
        1163             :       check_hist_coord_r8 = .false.
        1164             :     end if
        1165           0 :   end function check_hist_coord_r8
        1166             :
        1167           0 :   logical function check_hist_coord_r81(defined, input)
        1168             :
        1169             :     ! Input variables
        1170             :     real(r8),             pointer            :: defined(:)
        1171             :     real(r8), intent(in),          optional  :: input(:)
        1172             :
        1173             :     ! Local variables
        1174             :     integer                                  :: i
        1175             :
        1176           0 :     if (.not. associated(defined)) then
        1177             :       ! In this case, we assume the current value is undefined so any input OK
        1178             :       check_hist_coord_r81 = .true.
        1179           0 :     else if (present(input)) then
        1180             :       ! We have to match definitions
        1181           0 :       check_hist_coord_r81 = (size(input) == size(defined))
        1182             :     else
        1183             :       ! Not sure here. We have a value and are redefining without one?
        1184             :       check_hist_coord_r81 = .false.
        1185             :     end if
        1186           0 :     if (check_hist_coord_r81 .and. associated(defined)) then
        1187             :       ! Need to check the values
        1188           0 :       do i = 1, size(defined)
        1189           0 :         if (defined(i) /= input(i)) then
        1190             :           check_hist_coord_r81 = .false.
        1191             :           exit
        1192             :         end if
        1193             :       end do
        1194             :     end if
        1195           0 :   end function check_hist_coord_r81
        1196             :
        1197           0 :   logical function check_hist_coord_r82(defined, input)
        1198             :
        1199             :     ! Input variables
        1200             :     real(r8),             pointer            :: defined(:,:)
        1201             :     real(r8), intent(in),          optional  :: input(:,:)
        1202             :
        1203             :     ! Local variables
        1204             :     integer                                  :: i, j
        1205             :
        1206           0 :     if (.not. associated(defined)) then
        1207             :       ! In this case, we assume the current value is undefined so any input OK
        1208             :       check_hist_coord_r82 = .true.
        1209           0 :     else if (present(input)) then
        1210             :       ! We have to match definitions
        1211             :       check_hist_coord_r82 = ((size(input, 1) == size(defined, 1)) .and.    &
        1212           0 :                               (size(input, 2) == size(defined, 2)))
        1213             :     else
        1214             :       ! Not sure here. We have a value and are redefining without one?
        1215             :       check_hist_coord_r82 = .false.
        1216             :     end if
        1217           0 :     if (check_hist_coord_r82 .and. associated(defined)) then
        1218             :       ! Need to check the values
        1219           0 :       do j = 1, size(defined, 2)
        1220           0 :         do i = 1, size(defined, 1)
        1221           0 :           if (defined(i, j) /= input(i, j)) then
        1222             :             check_hist_coord_r82 = .false.
        1223             :             exit
        1224             :           end if
        1225             :         end do
        1226             :       end do
        1227             :     end if
        1228           0 :   end function check_hist_coord_r82
        1229             :
        1230           0 :   logical function check_hist_coord_ft(defined, input)
        1231             :
        1232             :     ! Input variables
        1233             :     type(formula_terms_t), intent(in)           :: defined
        1234             :     type(formula_terms_t), intent(in), optional :: input
        1235             :
        1236             :     ! We will assume that if formula_terms has been defined, a_name has a value
        1237           0 :     if (len_trim(defined%a_name) == 0) then
        1238             :       ! In this case, we assume the current value is undefined so any input OK
        1239             :       check_hist_coord_ft = .true.
        1240           0 :     else if (present(input)) then
        1241             :       ! We have to match definitions
        1242             :       ! Need to check the values
        1243             :       check_hist_coord_ft =                                                   &
        1244             :            check_hist_coord(defined%a_name,       input%a_name)         .and. &
        1245             :            check_hist_coord(defined%a_long_name,  input%a_long_name)    .and. &
        1246             :            check_hist_coord(defined%a_values,     input%a_values)       .and. &
        1247             :            check_hist_coord(defined%b_name,       input%b_name)         .and. &
        1248             :            check_hist_coord(defined%b_long_name,  input%b_long_name)    .and. &
        1249             :            check_hist_coord(defined%b_values,     input%b_values)       .and. &
        1250             :            check_hist_coord(defined%p0_name,      input%p0_name)        .and. &
        1251             :            check_hist_coord(defined%p0_long_name, input%p0_long_name)   .and. &
        1252             :            check_hist_coord(defined%p0_units,     input%p0_units)       .and. &
        1253             :            check_hist_coord(defined%p0_value,     input%p0_value)       .and. &
        1254           0 :            check_hist_coord(defined%ps_name,      input%ps_name)
        1255             :     else
        1256             :       ! Not sure here. We have a value and are redefining without one?
        1257             :       check_hist_coord_ft = .false.
        1258             :     end if
        1259           0 :   end function check_hist_coord_ft
        1260             :
        1261             :   ! check_hist_coord: returns 0 if <name> is not registered as a hist coord
        1262             :   !                   returns i if <name> is registered with compatible values
        1263             :   !                   calls endrun if <name> is registered with incompatible values
        1264       18432 :   integer function check_hist_coord_all(name, vlen, long_name, units, bounds, &
        1265       18432 :        i_values, r_values, bounds_name, positive, standard_name, formula_terms)
        1266             :
        1267             :     ! Input variables
        1268             :     character(len=*),      intent(in)            :: name
        1269             :     integer,               intent(in)            :: vlen
        1270             :     character(len=*),      intent(in),  optional :: long_name
        1271             :     character(len=*),      intent(in),  optional :: units
        1272             :     character(len=*),      intent(in),  optional :: bounds_name
        1273             :     integer,               intent(in),  optional :: i_values(:)
        1274             :     real(r8),              intent(in),  optional :: r_values(:)
        1275             :     real(r8),              intent(in),  optional :: bounds(:,:)
        1276             :     character(len=*),      intent(in),  optional :: positive
        1277             :     character(len=*),      intent(in),  optional :: standard_name
        1278             :     type(formula_terms_t), intent(in),  optional :: formula_terms
        1279             :
        1280             :     ! Local variables
        1281             :     character(len=120)                           :: errormsg
        1282             :     integer                                      :: i
        1283             :
        1284       18432 :     i = get_hist_coord_index(trim(name))
        1285             :     ! If i > 0, this mdim has already been registered
        1286       18432 :     if (i > 0) then
        1287           0 :       check_hist_coord_all = i
        1288           0 :       if (.not. check_hist_coord(hist_coords(i)%dimsize, vlen)) then
        1289           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, '//trim(name)//' with incompatible size'
        1290           0 :         call endrun(errormsg)
        1291             :       end if
        1292           0 :       if (.not. check_hist_coord(hist_coords(i)%long_name, long_name)) then
        1293           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with a different long_name'
        1294           0 :         call endrun(errormsg)
        1295             :       end if
        1296           0 :       if (.not. check_hist_coord(hist_coords(i)%units, units)) then
        1297           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with different units'
        1298           0 :         call endrun(errormsg)
        1299             :       end if
        1300           0 :       if (.not. check_hist_coord(hist_coords(i)%bounds_name, bounds_name)) then
        1301           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with a different bounds_name'
        1302           0 :         call endrun(errormsg)
        1303             :       end if
        1304           0 :       if (.not. check_hist_coord(hist_coords(i)%standard_name, standard_name)) then
        1305           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with a different standard_name'
        1306           0 :         call endrun(errormsg)
        1307             :       end if
        1308           0 :       if (.not. check_hist_coord(hist_coords(i)%positive, positive)) then
        1309           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with a different value of positive'
        1310           0 :         call endrun(errormsg)
        1311             :       end if
        1312             :       ! Since the integer_dim defaults to .true., double check which to check
        1313           0 :       if ((.not. hist_coords(i)%integer_dim) .or.                             &
        1314             :            associated(hist_coords(i)%real_values)) then
        1315           0 :         if (.not. check_hist_coord(hist_coords(i)%real_values, r_values)) then
        1316           0 :           write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with different values'
        1317           0 :           call endrun(errormsg)
        1318           0 :         else if (present(i_values)) then
        1319           0 :           write(errormsg, *) 'ERROR: Attempt to register integer values for real dimension'
        1320           0 :           call endrun(errormsg)
        1321             :         end if
        1322             :       else
        1323           0 :         if (.not. check_hist_coord(hist_coords(i)%integer_values, i_values)) then
        1324           0 :           write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with different values'
        1325           0 :           call endrun(errormsg)
        1326           0 :         else if (present(i_values) .and. present(r_values)) then
        1327           0 :           write(errormsg, *) 'ERROR: Attempt to register real values for integer dimension'
        1328           0 :           call endrun(errormsg)
        1329             :         end if
        1330             :       end if
        1331           0 :       if (.not. check_hist_coord(hist_coords(i)%bounds, bounds)) then
        1332           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with different bounds'
        1333           0 :         call endrun(errormsg)
        1334             :       end if
        1335           0 :       if (.not. check_hist_coord(hist_coords(i)%formula_terms, formula_terms)) then
        1336           0 :         write(errormsg, *) 'ERROR: Attempt to register dimension, ',trim(name),' with different formula_terms'
        1337           0 :         call endrun(errormsg)
        1338             :       end if
        1339             :     else
        1340             :       check_hist_coord_all = 0
        1341             :     end if
        1342       18432 :   end function check_hist_coord_all
        1343             :
        1344       12288 :   subroutine add_hist_coord_regonly(name, index)
        1345             :
        1346             :     ! Input variable
        1347             :     character(len=*),  intent(in)    :: name
        1348             :     integer, optional, intent(out)   :: index
        1349             :
        1350             :     ! Local variables
        1351             :     character(len=120)               :: errormsg
        1352             :     integer                          :: i
        1353             :
        1354       12288 :     if ((trim(name) == trim(horiz_only)) .or. (len_trim(name) == 0)) then
        1355           0 :       call endrun('ADD_HIST_COORD: '//trim(name)//' is not a valid coordinate name')
        1356             :     end if
        1357       12288 :     i = get_hist_coord_index(trim(name))
        1358             :     ! If i > 0, this mdim has already been registered
        1359       12288 :     if (i <= 0) then
        1360       12288 :       registeredmdims = registeredmdims + 1
        1361       12288 :       if (registeredmdims > maxmdims) then
        1362           0 :         call endrun('Too many dimensions in add_hist_coord.')
        1363             :       end if
        1364       12288 :       if (len_trim(name) > max_hcoordname_len) then
        1365           0 :         write(errormsg,'(a,i3,a)') 'History coord name exceeds the ',         &
        1366           0 :              max_hcoordname_len, ' character length limit'
        1367           0 :         call endrun(errormsg)
        1368             :       end if
        1369       12288 :       hist_coords(registeredmdims)%name = trim(name)
        1370       12288 :       hist_coords(registeredmdims)%dimsize = 0
        1371       12288 :       hist_coords(registeredmdims)%long_name = ''
        1372       12288 :       hist_coords(registeredmdims)%units = ''
        1373       12288 :       hist_coords(registeredmdims)%integer_dim = .true.
        1374       12288 :       hist_coords(registeredmdims)%positive = ''
        1375       12288 :       hist_coords(registeredmdims)%standard_name = ''
        1376       12288 :       if (present(index)) then
        1377       12288 :         index = registeredmdims
        1378             :       end if
        1379             :     else
        1380           0 :       if (present(index)) then
        1381           0 :         index = i
        1382             :       end if
        1383             :     end if
        1384             :
        1385       12288 :   end subroutine add_hist_coord_regonly
        1386             :
        1387           0 :   subroutine add_hist_coord_int(name, vlen, long_name, units, values,         &
        1388             :        positive, standard_name, dimname)
        1389             :
        1390             :     ! Input variables
        1391             :     character(len=*), intent(in)                    :: name
        1392             :     integer,          intent(in)                    :: vlen
        1393             :     character(len=*), intent(in)                    :: long_name
        1394             :     character(len=*), intent(in),          optional :: units
        1395             :     integer,          intent(in),  target, optional :: values(:)
        1396             :     character(len=*), intent(in),          optional :: positive
        1397             :     character(len=*), intent(in),          optional :: standard_name
        1398             :     character(len=*), intent(in),          optional :: dimname
        1399             :
        1400             :     ! Local variables
        1401             :     integer                                         :: i
        1402             :
        1403             :     ! First, check to see if it is OK to add this coord
        1404             :     i = check_hist_coord(name, vlen=vlen, long_name=long_name, units=units,   &
        1405           0 :          i_values=values, positive=positive, standard_name=standard_name)
        1406             :     ! Register the name if necessary
        1407           0 :     if (i == 0) then
        1408           0 :        call add_hist_coord(trim(name), i)
        1409           0 :        if(masterproc) then
        1410           0 :           write(iulog, '(3a,i0,a,i0)') 'Registering hist coord: ', trim(name),  &
        1411           0 :                '(', i, ') with length: ', vlen
        1412             :        end if
        1413             :     end if
        1414             :
        1415             :     ! Set the coord's values
        1416           0 :     hist_coords(i)%dimsize = vlen
        1417           0 :     if (len_trim(long_name) > max_chars) then
        1418           0 :        if(masterproc) then
        1419           0 :           write(iulog,*) 'WARNING: long_name for ',trim(name),' too long'
        1420             :        end if
        1421             :     end if
        1422           0 :     hist_coords(i)%long_name = trim(long_name)
        1423           0 :     if (present(units)) then
        1424           0 :        hist_coords(i)%units = trim(units)
        1425             :     else
        1426           0 :        hist_coords(i)%units = ''
        1427             :     end if
        1428           0 :     hist_coords(i)%integer_dim = .true.
        1429           0 :     if (present(values)) then
        1430           0 :        hist_coords(i)%integer_values => values
        1431             :     endif
        1432           0 :     if (present(positive)) then
        1433           0 :        hist_coords(i)%positive = trim(positive)
        1434             :     end if
        1435           0 :     if (present(standard_name)) then
        1436           0 :        hist_coords(i)%standard_name = trim(standard_name)
        1437             :     end if
        1438           0 :     hist_coords(i)%vertical_coord = .false.
        1439           0 :     if (present(dimname)) then
        1440           0 :        hist_coords(i)%dimname = trim(dimname)
        1441             :     else
        1442           0 :        hist_coords(i)%dimname = ''
        1443             :     end if
        1444             :
        1445           0 :   end subroutine add_hist_coord_int
        1446             :
        1447       12288 :   subroutine add_hist_coord_r8(name, vlen, long_name, units, values,         &
        1448       12288 :        bounds_name, bounds, positive, standard_name, vertical_coord, dimname)
        1449             :
        1450             :     ! Input variables
        1451             :     character(len=*),      intent(in)                    :: name
        1452             :     integer,               intent(in)                    :: vlen
        1453             :     character(len=*),      intent(in)                    :: long_name
        1454             :     character(len=*),      intent(in)                    :: units
        1455             :     real(r8),              intent(in), target            :: values(:)
        1456             :     character(len=*),      intent(in),          optional :: bounds_name
        1457             :     real(r8),              intent(in), target,  optional :: bounds(:,:)
        1458             :     character(len=*),      intent(in),          optional :: positive
        1459             :     character(len=*),      intent(in),          optional :: standard_name
        1460             :     logical,               intent(in),          optional :: vertical_coord
        1461             :     character(len=*),      intent(in),          optional :: dimname
        1462             :
        1463             :     ! Local variables
        1464             :     character(len=120)                                   :: errormsg
        1465             :     integer                                              :: i
        1466             :
        1467             :     ! First, check to see if it is OK to add this coord
        1468             :     i = check_hist_coord(name, vlen=vlen, long_name=long_name, units=units,   &
        1469             :          r_values=values, positive=positive, standard_name=standard_name,     &
        1470       55296 :          bounds_name=bounds_name, bounds=bounds)
        1471             :     ! Register the name if necessary
        1472       12288 :     if (i == 0) then
        1473       12288 :        call add_hist_coord(trim(name), i)
        1474       12288 :        if(masterproc) then
        1475          16 :           write(iulog, '(3a,i0,a,i0)') 'Registering hist coord: ', trim(name),  &
        1476          32 :                '(', i, ') with length: ', vlen
        1477             :        end if
        1478             :     end if
        1479             :
        1480             :     ! Set the coord's size
        1481       12288 :     hist_coords(i)%dimsize = vlen
        1482       12288 :     if (len_trim(long_name) > max_chars) then
        1483           0 :        if(masterproc) then
        1484           0 :           write(iulog,*) 'WARNING: long_name for ',trim(name),' too long'
        1485             :        end if
        1486             :     end if
        1487       12288 :     hist_coords(i)%long_name = trim(long_name)
        1488       12288 :     if (len_trim(units) > 0) then
        1489       12288 :        hist_coords(i)%units = trim(units)
        1490             :     else
        1491           0 :        hist_coords(i)%units = '1'
        1492             :     end if
        1493       12288 :     hist_coords(i)%integer_dim = .false.
        1494       12288 :     hist_coords(i)%real_values => values
        1495       12288 :     if (present(positive)) then
        1496        6144 :        hist_coords(i)%positive = trim(positive)
        1497             :     end if
        1498       12288 :     if (present(standard_name)) then
        1499           0 :        hist_coords(i)%standard_name = trim(standard_name)
        1500             :     end if
        1501       12288 :     if (present(bounds_name)) then
        1502           0 :        hist_coords(i)%bounds_name = trim(bounds_name)
        1503           0 :        if (.not. present(bounds)) then
        1504           0 :           write(errormsg,*) 'bounds must be present for ',trim(bounds_name)
        1505           0 :           call endrun(errormsg)
        1506             :        end if
        1507           0 :        hist_coords(i)%bounds => bounds
        1508       12288 :     else if (present(bounds)) then
        1509           0 :        write(errormsg,*) 'bounds_name must be present for bounds values'
        1510           0 :        call endrun(errormsg)
        1511             :     else
        1512       12288 :        hist_coords(i)%bounds_name = ''
        1513             :     end if
        1514       12288 :     if (present(vertical_coord)) then
        1515        6144 :        hist_coords(i)%vertical_coord = vertical_coord
        1516             :     else
        1517        6144 :        hist_coords(i)%vertical_coord = .false.
        1518             :     end if
        1519       12288 :     if (present(dimname)) then
        1520        6144 :        hist_coords(i)%dimname = trim(dimname)
        1521             :     else
        1522        6144 :        hist_coords(i)%dimname = ''
        1523             :     end if
        1524             :
        1525       12288 :   end subroutine add_hist_coord_r8
        1526             :
        1527        6144 :   subroutine add_vert_coord(name, vlen, long_name, units, values,            &
        1528             :        positive, standard_name, formula_terms)
        1529             :
        1530             :     ! Input variables
        1531             :     character(len=*),      intent(in)                    :: name
        1532             :     integer,               intent(in)                    :: vlen
        1533             :     character(len=*),      intent(in)                    :: long_name
        1534             :     character(len=*),      intent(in)                    :: units
        1535             :     real(r8),              intent(in), target            :: values(:)
        1536             :     character(len=*),      intent(in),          optional :: positive
        1537             :     character(len=*),      intent(in),          optional :: standard_name
        1538             :     type(formula_terms_t), intent(in),          optional :: formula_terms
        1539             :
        1540             :     ! Local variable
        1541             :     integer                                              :: i
        1542             :
        1543             :     ! First, check to see if it is OK to add this coord
        1544             :     i = check_hist_coord(name, vlen=vlen, long_name=long_name, units=units,   &
        1545             :          r_values=values, positive=positive, standard_name=standard_name,     &
        1546       12288 :          formula_terms=formula_terms)
        1547             :     ! Register the name and hist_coord values if necessary
        1548        6144 :     if (i == 0) then
        1549             :       call add_hist_coord(trim(name), vlen, long_name, units, values,         &
        1550             :            positive=positive, standard_name=standard_name,                    &
        1551        6144 :            vertical_coord=.true.)
        1552        6144 :       i = get_hist_coord_index(trim(name))
        1553        6144 :       if(masterproc) then
        1554           8 :          write(iulog, '(3a,i0,a,i0)') 'Registering hist coord: ', trim(name),   &
        1555          16 :               '(', i, ') with length: ', vlen
        1556             :       end if
        1557             :     end if
        1558             :
        1559        6144 :     if (present(formula_terms)) then
        1560           0 :       hist_coords(i)%formula_terms = formula_terms
        1561             :     end if
        1562             :
        1563        6144 :   end subroutine add_vert_coord
        1564             :
        1565     1966080 :   subroutine write_hist_coord_attr(File, mdimind, boundsdim, dimonly, mdimid)
        1566             :     use pio, only: file_desc_t, var_desc_t, pio_put_att, pio_noerr,           &
        1567             :                    pio_int, pio_double, pio_inq_varid, pio_def_var
        1568             :     use cam_pio_utils, only: cam_pio_def_dim, cam_pio_def_var
        1569             :
        1570             :     ! Input variables
        1571             :     type(file_desc_t), intent(inout) :: File           ! PIO file Handle
        1572             :     integer,           intent(in)    :: mdimind        ! Internal dim index
        1573             :     integer,           intent(in)    :: boundsdim      ! Bounds dimension ID
        1574             :     logical,           intent(in)    :: dimonly        ! No def_var if .true.
        1575             :     integer, optional, intent(out)   :: mdimid
        1576             :
        1577             :     ! Local variables
        1578             :     integer                          :: dimid          ! PIO dimension ID
        1579             :     type(var_desc_t)                 :: vardesc        ! PIO variable descriptor
        1580             :     character(len=120)               :: errormsg
        1581             :     character(len=max_chars)         :: formula_terms  ! Constructed string
        1582             :     integer                          :: ierr
        1583             :     integer                          :: dtype
        1584             :     logical                          :: defvar         ! True if var exists
        1585             :
        1586             :     ! Create or check dimension for this coordinate
        1587     1966080 :     if (len_trim(hist_coords(mdimind)%dimname) > 0) then
        1588             :        ! Dim can already exist if different from coord name
        1589             :        call cam_pio_def_dim(File, trim(hist_coords(mdimind)%dimname),         &
        1590             :             hist_coords(mdimind)%dimsize, dimid,                              &
        1591             :             existOK=(trim(hist_coords(mdimind)%dimname) /=                    &
        1592      983040 :                      trim(hist_coords(mdimind)%name)))
        1593             :     else
        1594             :        ! The dimension has the same name as the coord -- must be new dim
        1595             :        call cam_pio_def_dim(File, trim(hist_coords(mdimind)%name),            &
        1596      983040 :             hist_coords(mdimind)%dimsize, dimid, existOK=.false.)
        1597             :     end if
        1598             :     ! If the caller wants to know the NetCDF dimension ID, set it here
        1599     1966080 :     if (present(mdimid)) then
        1600     1966080 :       mdimid = dimid
        1601             :     end if
        1602     1966080 :     if (.not. dimonly) then
        1603             :       ! Time to define the variable (only if there are values)
        1604     1966080 :       if (hist_coords(mdimind)%integer_dim) then
        1605           0 :         dtype = pio_int
        1606           0 :         defvar = associated(hist_coords(mdimind)%integer_values)
        1607             :       else
        1608     1966080 :         dtype = pio_double
        1609     1966080 :         defvar = associated(hist_coords(mdimind)%real_values)
        1610             :       end if
        1611     1966080 :       if (defvar) then
        1612             :         call cam_pio_def_var(File, trim(hist_coords(mdimind)%name), dtype,    &
        1613     3932160 :              (/dimid/), vardesc, existOK=.false.)
        1614             :         ! long_name
        1615     1966080 :         ierr=pio_put_att(File, vardesc, 'long_name', trim(hist_coords(mdimind)%long_name))
        1616     1966080 :         call cam_pio_handle_error(ierr, 'Error writing "long_name" attr in write_hist_coord_attr')
        1617             :         ! units
        1618     1966080 :         if(len_trim(hist_coords(mdimind)%units) > 0) then
        1619             :           ierr=pio_put_att(File, vardesc, 'units',                              &
        1620     1966080 :                trim(hist_coords(mdimind)%units))
        1621     1966080 :           call cam_pio_handle_error(ierr, 'Error writing "units" attr in write_hist_coord_attr')
        1622             :         end if
        1623             :         ! positive
        1624     1966080 :         if(len_trim(hist_coords(mdimind)%positive) > 0) then
        1625             :           ierr=pio_put_att(File, vardesc, 'positive',                           &
        1626      983040 :                trim(hist_coords(mdimind)%positive))
        1627      983040 :           call cam_pio_handle_error(ierr, 'Error writing "positive" attr in write_hist_coord_attr')
        1628             :         end if
        1629             :         ! standard_name
        1630     1966080 :         if(len_trim(hist_coords(mdimind)%standard_name) > 0) then
        1631             :           ierr=pio_put_att(File, vardesc, 'standard_name',                      &
        1632           0 :                trim(hist_coords(mdimind)%standard_name))
        1633           0 :           call cam_pio_handle_error(ierr, 'Error writing "standard_name" attr in write_hist_coord_attr')
        1634             :         end if
        1635             :         ! formula_terms
        1636     1966080 :         if(len_trim(hist_coords(mdimind)%formula_terms%a_name) > 0) then
        1637             :           write(formula_terms, '("a: ",a," b: ",a," p0: ",a," ps: ",a)')        &
        1638           0 :                trim(hist_coords(mdimind)%formula_terms%a_name),                 &
        1639           0 :                trim(hist_coords(mdimind)%formula_terms%b_name),                 &
        1640           0 :                trim(hist_coords(mdimind)%formula_terms%p0_name),                &
        1641           0 :                trim(hist_coords(mdimind)%formula_terms%ps_name)
        1642           0 :           ierr=pio_put_att(File, vardesc, 'formula_terms', trim(formula_terms))
        1643           0 :           call cam_pio_handle_error(ierr, 'Error writing "formula_terms" attr in write_hist_coord_attr')
        1644             :         end if
        1645             :         ! bounds
        1646     1966080 :         if (associated(hist_coords(mdimind)%bounds)) then
        1647             :           ! Write name of the bounds variable
        1648           0 :           ierr=pio_put_att(File, vardesc, 'bounds', trim(hist_coords(mdimind)%bounds_name))
        1649           0 :           call cam_pio_handle_error(ierr, 'Error writing "bounds" attr in write_hist_coord_attr')
        1650             :         end if
        1651             :       end if
        1652             :
        1653             :       ! Now, we need to define and populate the associated bounds variable
        1654             :       ! NB: Reusing vardesc, no longer assocated with main variable
        1655     1966080 :       if (associated(hist_coords(mdimind)%bounds)) then
        1656           0 :         if (size(hist_coords(mdimind)%bounds,2) /= hist_coords(mdimind)%dimsize) then
        1657             :           ! If anyone hits this check, add a new dimension for this case
        1658           0 :           write(errormsg, *) 'The bounds variable, ',                         &
        1659           0 :                trim(hist_coords(mdimind)%bounds_name),                        &
        1660           0 :                ', needs to have dimension (2,', hist_coords(mdimind)%dimsize
        1661           0 :           call endrun(errormsg)
        1662             :         end if
        1663             :         call cam_pio_def_var(File, trim(hist_coords(mdimind)%bounds_name),    &
        1664           0 :              pio_double, (/boundsdim,dimid/), vardesc, existOK=.false.)
        1665             :       end if
        1666             :
        1667             :       ! See if we have formula_terms variables to define
        1668             :       ! Define the "a" variable name
        1669             :       ! NB: Reusing vardesc, no longer assocated with previous variables
        1670     1966080 :       if (associated(hist_coords(mdimind)%formula_terms%a_values)) then
        1671           0 :         if (size(hist_coords(mdimind)%formula_terms%a_values) /= hist_coords(mdimind)%dimsize) then
        1672           0 :           write(errormsg, *) 'The forumla_terms variable, ',                  &
        1673           0 :                trim(hist_coords(mdimind)%formula_terms%a_name),               &
        1674           0 :                ', needs to have dimension', hist_coords(mdimind)%dimsize
        1675           0 :           call endrun(errormsg)
        1676             :         end if
        1677             :         call cam_pio_def_var(File, trim(hist_coords(mdimind)%formula_terms%a_name), &
        1678           0 :              pio_double, (/dimid/), vardesc, existOK=.false.)
        1679           0 :         ierr = pio_put_att(File, vardesc, 'long_name', trim(hist_coords(mdimind)%formula_terms%a_long_name))
        1680           0 :         call cam_pio_handle_error(ierr, 'Error writing "long_name" attr for "a" formula_term in write_hist_coord_attr')
        1681             :       end if
        1682             :       ! Define the "b" variable name
        1683             :       ! NB: Reusing vardesc, no longer assocated with previous variables
        1684     1966080 :       if (associated(hist_coords(mdimind)%formula_terms%b_values)) then
        1685           0 :         if (size(hist_coords(mdimind)%formula_terms%b_values) /= hist_coords(mdimind)%dimsize) then
        1686           0 :           write(errormsg, *) 'The forumla_terms variable, ',                  &
        1687           0 :                trim(hist_coords(mdimind)%formula_terms%b_name),               &
        1688           0 :                ', needs to have dimension', hist_coords(mdimind)%dimsize
        1689           0 :           call endrun(errormsg)
        1690             :         end if
        1691             :         call cam_pio_def_var(File, trim(hist_coords(mdimind)%formula_terms%b_name), &
        1692           0 :              pio_double, (/dimid/), vardesc, existOK=.false.)
        1693           0 :         ierr = pio_put_att(File, vardesc, 'long_name', trim(hist_coords(mdimind)%formula_terms%b_long_name))
        1694           0 :         call cam_pio_handle_error(ierr, 'Error writing "long_name" attr for "b" formula_term in write_hist_coord_attr')
        1695             :       end if
        1696             :       ! Maybe define the p0 variable (this may be defined already which is OK)
        1697             :       ! NB: Reusing vardesc, no longer assocated with previous variables
        1698     1966080 :       if (hist_coords(mdimind)%formula_terms%p0_value /= fillvalue) then
        1699           0 :         ierr = pio_inq_varid(File, trim(hist_coords(mdimind)%formula_terms%p0_name), vardesc)
        1700           0 :         if (ierr /= PIO_NOERR) then
        1701             :           ierr = pio_def_var(File, trim(hist_coords(mdimind)%formula_terms%p0_name), &
        1702           0 :                pio_double, vardesc)
        1703           0 :           call cam_pio_handle_error(ierr, 'Unable to define "p0" formula_terms variable in write_hist_coord_attr')
        1704           0 :           ierr = pio_put_att(File, vardesc, 'long_name', trim(hist_coords(mdimind)%formula_terms%p0_long_name))
        1705           0 :           call cam_pio_handle_error(ierr, 'Error writing "long_name" attr for "p0" formula_term in write_hist_coord_attr')
        1706           0 :           ierr = pio_put_att(File, vardesc, 'units', trim(hist_coords(mdimind)%formula_terms%p0_units))
        1707           0 :           call cam_pio_handle_error(ierr, 'Error writing "units" attr for "p0" formula_term in write_hist_coord_attr')
        1708             :         end if
        1709             :       end if
        1710             :       ! PS is not our responsibility
        1711             :     end if ! (.not. dimonly)
        1712             :
        1713     1966080 :   end subroutine write_hist_coord_attr
        1714             :
        1715             :   !---------------------------------------------------------------------------
        1716             :   !
        1717             :   !  write_hist_coord_attrs
        1718             :   !
        1719             :   !  Write the dimension and coordinate attributes for the defined history
        1720             :   !  coordinates.
        1721             :   !
        1722             :   !---------------------------------------------------------------------------
        1723      245760 :   subroutine write_hist_coord_attrs(File, boundsdim, mdimids, writemdims_in)
        1724     1966080 :     use pio, only: file_desc_t, var_desc_t, pio_put_att,         &
        1725             :                    pio_bcast_error, pio_internal_error, pio_seterrorhandling, &
        1726             :                    pio_char
        1727             :     use cam_pio_utils, only: cam_pio_def_dim, cam_pio_def_var
        1728             :
        1729             :     ! Input variables
        1730             :     type(file_desc_t), intent(inout) :: File           ! PIO file Handle
        1731             :     integer,           intent(in)    :: boundsdim      ! Bounds dimension ID
        1732             :     integer, optional, allocatable, intent(out)   :: mdimids(:) ! NetCDF dim IDs
        1733             :     logical, optional, intent(in)    :: writemdims_in  ! Write mdim variable
        1734             :
        1735             :     ! Local variables
        1736             :     integer                          :: i
        1737             :     integer                          :: ierr
        1738             :     integer                          :: dimids(2)      ! PIO dimension IDs
        1739             :     logical                          :: writemdims     ! Define an mdim variable
        1740             :     type(var_desc_t)                 :: vardesc        ! PIO variable descriptor
        1741             :
        1742      245760 :     if (present(mdimids)) then
        1743      737280 :       allocate(mdimids(registeredmdims))
        1744             :     end if
        1745             :
        1746             :     ! We will handle errors for this routine
        1747      245760 :     call pio_seterrorhandling(File, PIO_BCAST_ERROR)
        1748             :
        1749      245760 :     if (present(writemdims_in)) then
        1750      245760 :       writemdims = writemdims_in
        1751             :     else
        1752           0 :       writemdims = .false.
        1753             :     end if
        1754             :
        1755             :     ! NB: Currently, writemdims is for restart and we don't need to write
        1756             :     ! these out in a history-restart file. This could change in the future.
        1757             :     ! which would require a change to the function of the fourth argument
        1758             :     ! Fill in the attribute information for each mdim
        1759     2211840 :     do i = 1, registeredmdims
        1760     2211840 :       if (present(mdimids)) then
        1761     1966080 :         call write_hist_coord_attr(File, i, boundsdim, writemdims, mdimids(i))
        1762             :       else
        1763           0 :         call write_hist_coord_attr(File, i, boundsdim, writemdims)
        1764             :       end if
        1765             :     end do
        1766             :
        1767      245760 :     if (writemdims) then
        1768             :       call cam_pio_def_dim(File, 'mdimslen', max_hcoordname_len, dimids(1),   &
        1769           0 :            existOK=.true.)
        1770             :       call cam_pio_def_dim(File, 'num_mdims', registeredmdims, dimids(2),     &
        1771           0 :            existOK=.true.)
        1772             :       call cam_pio_def_var(File, mdim_var_name, pio_char, dimids, vardesc,    &
        1773           0 :            existOK=.false.)
        1774           0 :       ierr = pio_put_att(File, vardesc, 'long_name', 'mdim dimension names')
        1775           0 :       call cam_pio_handle_error(ierr, 'Error writing "long_name" attr for mdimnames in write_hist_coord_attrs')
        1776             :
        1777             :     end if
        1778             :
        1779             :     ! Back to I/O or die trying
        1780      245760 :     call pio_seterrorhandling(File, PIO_INTERNAL_ERROR)
        1781      245760 :   end subroutine write_hist_coord_attrs
        1782             :
        1783     1966080 :   subroutine write_hist_coord_var(File, mdimind)
        1784      245760 :     use pio, only: file_desc_t, var_desc_t, pio_put_var, pio_inq_varid
        1785             :
        1786             :     ! Input variables
        1787             :     type(file_desc_t), intent(inout) :: File           ! PIO file Handle
        1788             :     integer,           intent(in)    :: mdimind        ! Internal dim index
        1789             :
        1790             :     ! Local variables
        1791             :     type(var_desc_t)                 :: vardesc        ! PIO variable descriptor
        1792             :     integer                          :: ierr
        1793             :
        1794     1966080 :     if ((hist_coords(mdimind)%integer_dim .and.                               &
        1795     3932160 :          associated(hist_coords(mdimind)%integer_values)) .or.                &
        1796             :          ((.not. hist_coords(mdimind)%integer_dim) .and.                      &
        1797             :          associated(hist_coords(mdimind)%real_values))) then
        1798             :       ! Check to make sure the variable already exists in the file
        1799     1966080 :       ierr = pio_inq_varid(File, trim(hist_coords(mdimind)%name), vardesc)
        1800     1966080 :       call cam_pio_handle_error(ierr, 'Error writing values for nonexistent dimension variable write_hist_coord_var')
        1801             :       ! Write out the values for this dimension variable
        1802     1966080 :       if (hist_coords(mdimind)%integer_dim) then
        1803           0 :         ierr = pio_put_var(File, vardesc, hist_coords(mdimind)%integer_values)
        1804             :       else
        1805     1966080 :         ierr = pio_put_var(File, vardesc, hist_coords(mdimind)%real_values)
        1806             :       end if
        1807     1966080 :       call cam_pio_handle_error(ierr, 'Error writing variable values in write_hist_coord_var')
        1808             :     end if
        1809             :
        1810             :     ! Now, we need to possibly write values for the associated bounds variable
        1811     1966080 :     if (associated(hist_coords(mdimind)%bounds)) then
        1812             :       ! Check to make sure the variable already exists in the file
        1813             :       ! NB: Reusing vardesc, no longer assocated with previous variables
        1814           0 :       ierr = pio_inq_varid(File, trim(hist_coords(mdimind)%bounds_name), vardesc)
        1815           0 :       call cam_pio_handle_error(ierr, 'Error writing values for nonexistent bounds variable write_hist_coord_var')
        1816             :     ! Write out the values for this bounds variable
        1817           0 :       ierr = pio_put_var(File, vardesc, hist_coords(mdimind)%bounds)
        1818           0 :       call cam_pio_handle_error(ierr, 'Error writing bounds values in write_hist_coord_var')
        1819             :     end if
        1820             :
        1821             :     ! Write values for the "a" variable name
        1822     1966080 :     if (associated(hist_coords(mdimind)%formula_terms%a_values)) then
        1823             :       ! Check to make sure the variable already exists in the file
        1824             :       ! NB: Reusing vardesc, no longer assocated with previous variables
        1825           0 :       ierr = pio_inq_varid(File, trim(hist_coords(mdimind)%formula_terms%a_name), vardesc)
        1826           0 :       call cam_pio_handle_error(ierr, 'Error writing values for nonexistent "a" formula_terms variable write_hist_coord_var')
        1827             :     ! Write out the values for this "a" formula_terms variable
        1828           0 :       ierr = pio_put_var(File, vardesc, hist_coords(mdimind)%formula_terms%a_values)
        1829           0 :       call cam_pio_handle_error(ierr, 'Error writing "a" formula_terms values in write_hist_coord_var')
        1830             :     end if
        1831             :     ! Write values for the "b" variable name
        1832     1966080 :     if (associated(hist_coords(mdimind)%formula_terms%b_values)) then
        1833             :       ! Check to make sure the variable already exists in the file
        1834             :       ! NB: Reusing vardesc, no longer assocated with previous variables
        1835           0 :       ierr = pio_inq_varid(File, trim(hist_coords(mdimind)%formula_terms%b_name), vardesc)
        1836           0 :       call cam_pio_handle_error(ierr, 'Error writing values for nonexistent "b" formula_terms variable write_hist_coord_var')
        1837             :     ! Write out the values for this "b" formula_terms variable
        1838           0 :       ierr = pio_put_var(File, vardesc, hist_coords(mdimind)%formula_terms%b_values)
        1839           0 :       call cam_pio_handle_error(ierr, 'Error writing "b" formula_terms values in write_hist_coord_var')
        1840             :     end if
        1841             :     ! Write values for the "p0" variable name (this may be an overwrite, too bad)
        1842     1966080 :     if (hist_coords(mdimind)%formula_terms%p0_value /= fillvalue) then
        1843             :       ! Check to make sure the variable already exists in the file
        1844             :       ! NB: Reusing vardesc, no longer assocated with previous variables
        1845           0 :       ierr = pio_inq_varid(File, trim(hist_coords(mdimind)%formula_terms%p0_name), vardesc)
        1846           0 :       call cam_pio_handle_error(ierr, 'Error writing values for nonexistent "p0" formula_terms variable write_hist_coord_var')
        1847             :     ! Write out the values for this "p0" formula_terms variable
        1848           0 :       ierr = pio_put_var(File, vardesc, hist_coords(mdimind)%formula_terms%p0_value)
        1849           0 :       call cam_pio_handle_error(ierr, 'Error writing "p0" formula_terms values in write_hist_coord_var')
        1850             :     end if
        1851             :
        1852     1966080 :   end subroutine write_hist_coord_var
        1853             :
        1854      245760 :   subroutine write_hist_coord_vars(File, writemdims_in)
        1855             :    use pio, only: file_desc_t, var_desc_t, pio_put_var,                       &
        1856             :                   pio_bcast_error, pio_internal_error,                        &
        1857             :                   pio_seterrorhandling, pio_inq_varid
        1858             :
        1859             :     ! Input variables
        1860             :     type(file_desc_t), intent(inout) :: File           ! PIO file Handle
        1861             :     logical, optional, intent(in)    :: writemdims_in  ! Write mdim variable
        1862             :
        1863             :     ! Local variables
        1864             :     integer                          :: i
        1865             :     integer                          :: ierr
        1866             :     logical                          :: writemdims     ! Define an mdim variable
        1867             :     type(var_desc_t)                 :: vardesc        ! PIO variable descriptor
        1868      245760 :     character(len=max_hcoordname_len), allocatable :: mdimnames(:)
        1869             :
        1870             :     ! We will handle errors for this routine
        1871      245760 :     call pio_seterrorhandling(File, PIO_BCAST_ERROR)
        1872             :
        1873      245760 :     if (present(writemdims_in)) then
        1874      245760 :       writemdims = writemdims_in
        1875             :     else
        1876             :       writemdims = .false.
        1877             :     end if
        1878             :
        1879      245760 :     if (writemdims) then
        1880           0 :       allocate(mdimnames(registeredmdims))
        1881             :     end if
        1882             :
        1883             :     ! Write out the variable values for each mdim
        1884     2211840 :     do i = 1, registeredmdims
        1885     1966080 :       if (.not. writemdims) then
        1886             :         ! NB: Currently, writemdims is for restart and we don't need to write
        1887             :         ! these out in a history-restart file. This could change in the future
        1888             :         ! which is why it is a separate if block
        1889             :         ! Fill in the attribute information for each mdim
        1890     1966080 :         call write_hist_coord_var(File, i)
        1891             :       end if
        1892      245760 :       if (writemdims) then
        1893           0 :         mdimnames(i) = trim(hist_coords(i)%name)
        1894             :       end if
        1895             :     end do
        1896             :
        1897      245760 :     if (writemdims) then
        1898           0 :       ierr = pio_inq_varid(File, mdim_var_name, vardesc)
        1899           0 :       call cam_pio_handle_error(ierr, 'Error writing values for nonexistent mdimnames variable in write_hist_coord_vars')
        1900             :       ! Write out the values for mdim names
        1901           0 :       ierr = pio_put_var(File, vardesc, mdimnames)
        1902           0 :       call cam_pio_handle_error(ierr, 'Error writing values for mdimnames variable in write_hist_coord_vars')
        1903           0 :       deallocate(mdimnames)
        1904             :     end if
        1905             :
        1906             :     ! Back to I/O or die trying
        1907      245760 :     call pio_seterrorhandling(File, PIO_INTERNAL_ERROR)
        1908             :
        1909      245760 :   end subroutine write_hist_coord_vars
        1910             :
        1911      978432 :   subroutine lookup_hist_coord_indices(mdimnames, mdimindicies)
        1912             :     ! Dummy arguments
        1913             :     character(len=*), intent(in) :: mdimnames(:)
        1914             :     integer, intent(out) :: mdimindicies(:)
        1915             :
        1916             :     ! Local variables
        1917             :     integer :: i, j
        1918             :     integer :: cnt
        1919             :     character(len=120) :: errormsg
        1920             :     character(len=16) :: name
        1921             :
        1922             :
        1923      978432 :     cnt = size(mdimnames)
        1924     1592832 :     mdimindicies = -1
        1925             :
        1926             :
        1927     1592832 :     do j=1,cnt
        1928      614400 :       name = mdimnames(j)
        1929     6508032 :       do i = 1, registeredmdims
        1930     5529600 :         if(name .eq. hist_coords(i)%name) then
        1931      614400 :           mdimindicies(j)=i
        1932             :         end if
        1933             :       end do
        1934             :     end do
        1935     1592832 :     do j = 1, cnt
        1936     1592832 :       if(mdimindicies(j) < 0) then
        1937           0 :         do i = 1, registeredmdims
        1938           0 :           print *,__FILE__,__LINE__,i,hist_coords(i)%name
        1939             :         end do
        1940           0 :         write(errormsg,*) 'Name ',mdimnames(j),' is not a registered history coordinate'
        1941           0 :         call endrun(errormsg)
        1942             :       end if
        1943             :     end do
        1944             :
        1945      978432 :   end subroutine lookup_hist_coord_indices
        1946             :
        1947             :   ! Find the vertical dimension (if present) in dimnames and return its size
        1948             :   !    (which is the number of levels). Return -1 if not found
        1949             :   !    If dimnames is not present, search all of the registered history coords
        1950      978432 :   integer function hist_coord_find_levels(dimnames) result(levels)
        1951             :     ! Dummy argument
        1952             :     character(len=*), optional, intent(in) :: dimnames(:)
        1953             :
        1954             :     ! Local variables
        1955             :     integer i, index, dimcnt
        1956             :
        1957      978432 :     levels = -1  ! Error return value
        1958             :
        1959      978432 :     if (present(dimnames)) then
        1960      978432 :       dimcnt = size(dimnames)
        1961             :     else
        1962           0 :       dimcnt = registeredmdims
        1963             :     end if
        1964             :
        1965      978432 :     do i = 1, dimcnt
        1966      614400 :       if (present(dimnames)) then
        1967      614400 :         index = get_hist_coord_index(trim(dimnames(i)))
        1968      614400 :         if (i < 0) then
        1969           0 :           call endrun('hist_coord_find_levels: '//trim(dimnames(i))//' is not a registred history coordinate')
        1970             :         end if
        1971             :       else
        1972             :         index = i  ! Just cycle through all the registered mdims
        1973             :       end if
        1974             :
        1975      978432 :       if (hist_coords(index)%vertical_coord) then
        1976      614400 :         levels = hist_coords(index)%dimsize
        1977      614400 :         exit
        1978             :       end if
        1979             :     end do
        1980             :
        1981      978432 :   end function hist_coord_find_levels
        1982             :
        1983             :   !#######################################################################
        1984             :
        1985           0 :   character(len=max_hcoordname_len) function hist_dimension_name(size)
        1986             :   ! Given a specific size value, return the first registered dimension name which matches the size, if it exists
        1987             :   ! Otherwise the name returned is blank
        1988             :
        1989             :      integer, intent(in) :: size
        1990             :
        1991             :      integer :: i
        1992             :
        1993           0 :      hist_dimension_name = ''
        1994             :
        1995           0 :      do i=1, registeredmdims
        1996           0 :         if(size == hist_coords(i)%dimsize) then
        1997           0 :            hist_dimension_name = hist_coords(i)%name
        1998           0 :            exit
        1999             :         end if
        2000             :      end do
        2001             :
        2002           0 :   end function hist_dimension_name
        2003             :
        2004             :   !#######################################################################
        2005             :
        2006           0 :   subroutine hist_dimension_values_r8(name, rvalues, istart, istop, found)
        2007             :      ! Given the name of a dimension, return its (real) values in <rvalues>
        2008             :      ! If <istart> and <istop> are present, they are the beginning and ending
        2009             :      ! indices of the dimension values to return in <rvalues>. By default,
        2010             :      ! the entire array is copied.
        2011             :      ! If <found> is passed, return .true. if <name> is a defined dimension
        2012             :      !       with real values.
        2013             :
        2014             :      ! Dummy arguments
        2015             :      character(len=*),  intent(in)  :: name
        2016             :      real(r8),          intent(out) :: rvalues(:)
        2017             :      integer, optional, intent(in)  :: istart
        2018             :      integer, optional, intent(in)  :: istop
        2019             :      logical, optional, intent(out) :: found
        2020             :      ! Local variables
        2021             :      integer :: indx, jndx, rndx
        2022             :      integer :: ibeg
        2023             :      integer :: iend
        2024             :      logical :: dim_ok
        2025             :      real(r8), parameter :: unset_r8 = huge(1.0_r8)
        2026             :      character(len=*), parameter  :: subname = ': hist_dimension_values_r8'
        2027             :
        2028           0 :      dim_ok = .false.
        2029           0 :      rvalues(:) = unset_r8
        2030             :
        2031           0 :      do indx = 1, registeredmdims
        2032           0 :         if(trim(name) == trim(hist_coords(indx)%name)) then
        2033           0 :            dim_ok = associated(hist_coords(indx)%real_values)
        2034           0 :            if (dim_ok) then
        2035           0 :               if (present(istart)) then
        2036           0 :                  ibeg = istart
        2037           0 :                  if (ibeg < LBOUND(hist_coords(indx)%real_values, 1)) then
        2038           0 :                     call endrun(subname//": istart is outside the bounds")
        2039             :                  end if
        2040             :               else
        2041           0 :                  ibeg = LBOUND(hist_coords(indx)%real_values, 1)
        2042             :               end if
        2043           0 :               if (present(istop)) then
        2044           0 :                  iend = istop
        2045           0 :                  if (iend > UBOUND(hist_coords(indx)%real_values, 1)) then
        2046           0 :                     call endrun(subname//": istop is outside the bounds")
        2047             :                  end if
        2048             :               else
        2049           0 :                  iend = UBOUND(hist_coords(indx)%real_values, 1)
        2050             :               end if
        2051           0 :               if (SIZE(rvalues) < (iend - ibeg + 1)) then
        2052           0 :                  call endrun(subname//": rvalues too small")
        2053             :               end if
        2054           0 :               rndx = 1
        2055           0 :               do jndx = ibeg, iend
        2056           0 :                  rvalues(rndx) = hist_coords(indx)%real_values(jndx)
        2057           0 :                  rndx = rndx + 1
        2058             :               end do
        2059             :            end if
        2060           0 :            exit
        2061             :         end if
        2062             :      end do
        2063           0 :      if (present(found)) then
        2064           0 :         found = dim_ok
        2065             :      end if
        2066             :
        2067           0 :   end subroutine hist_dimension_values_r8
        2068             :
        2069             :   !#######################################################################
        2070             :
        2071           0 :   subroutine hist_dimension_values_int(name, ivalues, istart, istop, found)
        2072             :      ! Given the name of a dimension, return its (integer) values in <ivalues>
        2073             :      ! If <istart> and <istop> are present, they are the beginning and ending
        2074             :      ! indices of the dimension values to return in <ivalues>. By default,
        2075             :      ! the entire array is copied.
        2076             :      ! If <found> is passed, return .true. if <name> is a defined dimension
        2077             :      !       with integer values.
        2078             :
        2079             :      ! Dummy arguments
        2080             :      character(len=*),  intent(in)  :: name
        2081             :      integer,           intent(out) :: ivalues(:)
        2082             :      integer, optional, intent(in)  :: istart
        2083             :      integer, optional, intent(in)  :: istop
        2084             :      logical, optional, intent(out) :: found
        2085             :      ! Local variables
        2086             :      integer :: indx, jndx, rndx
        2087             :      integer :: ibeg
        2088             :      integer :: iend
        2089             :      logical :: dim_ok
        2090             :      integer, parameter  :: unset_i = huge(1)
        2091             :      character(len=*), parameter  :: subname = 'hist_dimension_values_int'
        2092             :
        2093           0 :      dim_ok = .false.
        2094           0 :      ivalues(:) = unset_i
        2095             :
        2096           0 :      do indx = 1, registeredmdims
        2097           0 :         if(trim(name) == trim(hist_coords(indx)%name)) then
        2098           0 :            dim_ok = associated(hist_coords(indx)%integer_values)
        2099           0 :            if (dim_ok) then
        2100           0 :               if (present(istart)) then
        2101           0 :                  ibeg = istart
        2102           0 :                  if (ibeg < LBOUND(hist_coords(indx)%integer_values, 1)) then
        2103           0 :                     call endrun(subname//": istart is outside the bounds")
        2104             :                  end if
        2105             :               else
        2106           0 :                  ibeg = LBOUND(hist_coords(indx)%integer_values, 1)
        2107             :               end if
        2108           0 :               if (present(istop)) then
        2109           0 :                  iend = istop
        2110           0 :                  if (iend > UBOUND(hist_coords(indx)%integer_values, 1)) then
        2111           0 :                     call endrun(subname//": istop is outside the bounds")
        2112             :                  end if
        2113             :               else
        2114           0 :                  iend = UBOUND(hist_coords(indx)%integer_values, 1)
        2115             :               end if
        2116           0 :               if (SIZE(ivalues) < (iend - ibeg + 1)) then
        2117           0 :                  call endrun(subname//": ivalues too small")
        2118             :               end if
        2119           0 :               rndx = 1
        2120           0 :               do jndx = ibeg, iend
        2121           0 :                  ivalues(rndx) = hist_coords(indx)%integer_values(jndx)
        2122           0 :                  rndx = rndx + 1
        2123             :               end do
        2124             :            end if
        2125           0 :            exit
        2126             :         end if
        2127             :      end do
        2128           0 :      if (present(found)) then
        2129           0 :         found = dim_ok
        2130             :      end if
        2131             :
        2132           0 :   end subroutine hist_dimension_values_int
        2133             :
        2134             :   !#######################################################################
        2135             :
        2136           0 : end module cam_history_support

Generated by: LCOV version 1.14