LCOV - code coverage report
Current view: | top level - chemistry/utils - tracer_data.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 571 | 1208 | 47.3 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 16 | 33 | 48.5 %

              Line data    Source code
           1             : module tracer_data
           2             : !-----------------------------------------------------------------------
           3             : ! module used to read (and interpolate) offline tracer data (sources and
           4             : ! mixing ratios)
           5             : ! Created by: Francis Vitt -- 2 May 2006
           6             : ! Modified by : Jim Edwards -- 10 March 2009
           7             : ! Modified by : Cheryl Craig and Chih-Chieh (Jack) Chen  -- February 2010
           8             : !-----------------------------------------------------------------------
           9             :
          10             :   use perf_mod,         only : t_startf, t_stopf
          11             :   use shr_kind_mod,     only : r8 => shr_kind_r8, shr_kind_cl
          12             :   use time_manager,     only : get_curr_date, get_step_size
          13             :   use spmd_utils,       only : masterproc
          14             :   use ppgrid,           only : pcols, pver, pverp, begchunk, endchunk
          15             :   use cam_abortutils,   only : endrun
          16             :   use cam_logfile,      only : iulog
          17             :
          18             :   use physics_buffer,   only : physics_buffer_desc, pbuf_get_field, pbuf_get_index
          19             :   use time_manager,     only : set_time_float_from_date, set_date_from_time_float
          20             :   use pio,              only : file_desc_t, var_desc_t, &
          21             :                                pio_seterrorhandling, pio_internal_error, pio_bcast_error, &
          22             :                                pio_char, pio_noerr, &
          23             :                                pio_inq_dimid, pio_inq_varid, &
          24             :                                pio_def_dim, pio_def_var, &
          25             :                                pio_put_att, pio_put_var, &
          26             :                                pio_get_var, pio_get_att, pio_nowrite, pio_inq_dimlen, &
          27             :                                pio_inq_vardimid, pio_inq_dimlen, pio_closefile, &
          28             :                                pio_inquire_variable
          29             :
          30             :   implicit none
          31             :
          32             :   private  ! all unless made public
          33             :   save
          34             :
          35             :   public :: trfld, input3d, input2d, trfile
          36             :   public :: trcdata_init
          37             :   public :: advance_trcdata
          38             :   public :: get_fld_data
          39             :   public :: get_fld_ndx
          40             :   public :: write_trc_restart
          41             :   public :: read_trc_restart
          42             :   public :: init_trc_restart
          43             :   public :: incr_filename
          44             :
          45             :
          46             :   ! !PUBLIC MEMBERS
          47             :
          48             :   type input3d
          49             :      real(r8), dimension(:,:,:), pointer :: data => null()
          50             :   endtype input3d
          51             :
          52             :   type input2d
          53             :      real(r8), dimension(:,:), pointer :: data => null()
          54             :   endtype input2d
          55             :
          56             :   type trfld
          57             :      real(r8), dimension(:,:,:), pointer :: data => null()
          58             :      type(input3d), dimension(4) :: input
          59             :      character(len=32) :: srcnam
          60             :      character(len=32) :: fldnam
          61             :      character(len=32) :: units
          62             :      type(var_desc_t) :: var_id
          63             :      integer :: coords(4) ! LATDIM | LONDIM | LEVDIM | TIMDIM
          64             :      integer :: order(4) ! LATDIM | LONDIM | LEVDIM | TIMDIM
          65             :      logical :: srf_fld = .false.
          66             :      integer :: pbuf_ndx = -1
          67             :   endtype trfld
          68             :
          69             :   type trfile
          70             :      type(input2d), dimension(4) :: ps_in
          71             :      character(len=shr_kind_cl) :: pathname = ' '
          72             :      character(len=shr_kind_cl) :: curr_filename = ' '
          73             :      character(len=shr_kind_cl) :: next_filename = ' '
          74             :      type(file_desc_t) :: curr_fileid
          75             :      type(file_desc_t) :: next_fileid
          76             :
          77             :      type(var_desc_t), pointer :: currfnameid => null() ! pio restart file var id
          78             :      type(var_desc_t), pointer :: nextfnameid => null() ! pio restart file var id
          79             :
          80             :      character(len=shr_kind_cl) :: filenames_list = ''
          81             :      real(r8) :: datatimem = -1.e36_r8     ! time of prv. values read in
          82             :      real(r8) :: datatimep = -1.e36_r8     ! time of nxt. values read in
          83             :      real(r8) :: datatimes(4)
          84             :      integer :: interp_recs
          85             :      real(r8), pointer, dimension(:) :: curr_data_times => null()
          86             :      real(r8), pointer, dimension(:) :: next_data_times => null()
          87             :      logical :: remove_trc_file = .false.  ! delete file when finished with it
          88             :      real(r8) :: offset_time
          89             :      integer :: cyc_ndx_beg
          90             :      integer :: cyc_ndx_end
          91             :      integer :: cyc_yr = 0
          92             :      real(r8) :: one_yr = 0
          93             :      real(r8) :: curr_mod_time ! model time - calendar day
          94             :      real(r8) :: next_mod_time ! model time - calendar day - next time step
          95             :      integer :: nlon = 0
          96             :      integer :: nlat = 0
          97             :      integer :: nlev = 0
          98             :      integer :: nilev = 0
          99             :      integer :: ps_coords(3) ! LATDIM | LONDIM | TIMDIM
         100             :      integer :: ps_order(3) ! LATDIM | LONDIM | TIMDIM
         101             :      real(r8), pointer, dimension(:) :: lons => null()
         102             :      real(r8), pointer, dimension(:) :: lats => null()
         103             :      real(r8), pointer, dimension(:) :: levs => null()
         104             :      real(r8), pointer, dimension(:) :: ilevs => null()
         105             :      real(r8), pointer, dimension(:) :: hyam => null()
         106             :      real(r8), pointer, dimension(:) :: hybm => null()
         107             :      real(r8), pointer, dimension(:) :: hyai => null()
         108             :      real(r8), pointer, dimension(:) :: hybi => null()
         109             :      real(r8), pointer, dimension(:,:) :: weight_x => null(), weight_y => null()
         110             :      integer, pointer, dimension(:) :: count_x => null(), count_y => null()
         111             :      integer, pointer, dimension(:,:) :: index_x => null(), index_y => null()
         112             :
         113             :      real(r8), pointer, dimension(:,:) :: weight0_x=>null(), weight0_y=>null()
         114             :      integer, pointer, dimension(:) :: count0_x=>null(), count0_y=>null()
         115             :      integer, pointer, dimension(:,:) :: index0_x=>null(), index0_y=>null()
         116             :      logical :: dist
         117             :
         118             :      real(r8)                        :: p0
         119             :      type(var_desc_t) :: ps_id
         120             :      logical,  allocatable, dimension(:) :: in_pbuf
         121             :      logical :: has_ps = .false.
         122             :      logical :: zonal_ave = .false.
         123             :      logical :: unstructured = .false.
         124             :      logical :: alt_data = .false.
         125             :      logical :: geop_alt = .false.
         126             :      logical :: cyclical = .false.
         127             :      logical :: cyclical_list = .false.
         128             :      logical :: weight_by_lat = .false.
         129             :      logical :: conserve_column = .false.
         130             :      logical :: fill_in_months = .false.
         131             :      logical :: fixed = .false.
         132             :      logical :: initialized = .false.
         133             :      logical :: top_bndry = .false.
         134             :      logical :: top_layer = .false.
         135             :      logical :: stepTime = .false.  ! Do not interpolate in time, but use stepwise times
         136             :   endtype trfile
         137             :
         138             :   integer, public, parameter :: MAXTRCRS = 100
         139             :
         140             :   integer, parameter :: LONDIM = 1
         141             :   integer, parameter :: LATDIM = 2
         142             :   integer, parameter :: LEVDIM = 3
         143             :   integer, parameter :: TIMDIM = 4
         144             :
         145             :   integer, parameter :: PS_TIMDIM = 3
         146             :
         147             :   integer, parameter :: ZA_LATDIM = 1
         148             :   integer, parameter :: ZA_LEVDIM = 2
         149             :   integer, parameter :: ZA_TIMDIM = 3
         150             :
         151             :   integer, parameter :: nm=1    ! array index for previous (minus) data
         152             :   integer, parameter :: np=2    ! array index for next (plus) data
         153             :
         154             :   integer :: plon, plat
         155             :
         156             :   integer,allocatable :: lon_global_grid_ndx(:,:)
         157             :   integer,allocatable :: lat_global_grid_ndx(:,:)
         158             :
         159             : contains
         160             :
         161             : !--------------------------------------------------------------------------
         162             : !--------------------------------------------------------------------------
         163        6144 :   subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
         164             :                            rmv_file, data_cycle_yr, data_fixed_ymd, data_fixed_tod, data_type )
         165             :
         166             :     use mo_constants,    only : d2r
         167             :     use dyn_grid,        only : get_dyn_grid_parm, get_horiz_grid_d
         168             :     use phys_grid,       only : get_rlat_all_p, get_rlon_all_p, get_ncols_p
         169             :     use dycore,          only : dycore_is
         170             :     use horizontal_interpolate, only : xy_interp_init
         171             :     use spmd_utils,       only: mpicom, mstrid=>masterprocid, mpi_real8, mpi_integer
         172             :
         173             :     implicit none
         174             :
         175             :     character(len=*),    intent(in)    :: specifier(:)
         176             :     character(len=*),    intent(in)    :: filename
         177             :     character(len=*),    intent(in)    :: filelist
         178             :     character(len=*),    intent(in)    :: datapath
         179             :     type(trfld), dimension(:), pointer :: flds
         180             :     type(trfile),        intent(inout) :: file
         181             :     logical,             intent(in)    :: rmv_file
         182             :     integer,             intent(in)    :: data_cycle_yr
         183             :     integer,             intent(in)    :: data_fixed_ymd
         184             :     integer,             intent(in)    :: data_fixed_tod
         185             :     character(len=*),    intent(in)    :: data_type
         186             :
         187             :     character(len=*), parameter :: sub = 'trcdata_init'
         188             :
         189             :     integer :: f, mxnflds, astat
         190             :     integer :: str_yr, str_mon, str_day
         191             :     integer :: lon_dimid, lat_dimid, lev_dimid, tim_dimid, old_dimid
         192             :     integer :: dimids(4), did
         193             :     type(var_desc_t) :: varid
         194             :     integer :: idx
         195             :     integer :: ierr
         196             :     integer :: errcode
         197             :     real(r8) :: start_time, time1, time2
         198             :     integer :: i1,i2,j1,j2
         199             :     integer :: nvardims, vardimids(4)
         200             :
         201             :     character(len=256) :: data_units
         202        6144 :     real(r8), allocatable :: lam(:), phi(:)
         203             :     real(r8):: rlats(pcols), rlons(pcols)
         204             :     integer :: lchnk, ncol, icol, i,j
         205             :     logical :: found
         206             :     integer :: aircraft_cnt
         207             :     integer :: err_handling
         208             :
         209        6144 :     call specify_fields( specifier, flds )
         210             :
         211        6144 :     file%datatimep=-1.e36_r8
         212        6144 :     file%datatimem=-1.e36_r8
         213             :
         214        6144 :     mxnflds = 0
         215        6144 :     if (associated(flds)) mxnflds = size( flds )
         216             :
         217        6144 :     if (mxnflds < 1) return
         218             :
         219        6144 :     file%remove_trc_file = rmv_file
         220        6144 :     file%pathname = trim(datapath)
         221        6144 :     file%filenames_list = trim(filelist)
         222             :
         223        6144 :     file%fill_in_months = .false.
         224        6144 :     file%cyclical = .false.
         225        6144 :     file%cyclical_list = .false.
         226        6144 :     file%dist = .false.
         227             :
         228           0 :     select case ( data_type )
         229             :     case( 'FIXED' )
         230           0 :        file%fixed = .true.
         231             :     case( 'INTERP_MISSING_MONTHS' )
         232           0 :        file%fill_in_months = .true.
         233             :     case( 'CYCLICAL' )
         234        3072 :        file%cyclical = .true.
         235        3072 :        file%cyc_yr = data_cycle_yr
         236             :     case( 'CYCLICAL_LIST' )
         237           0 :        file%cyclical_list = .true.
         238           0 :        file%cyc_yr = data_cycle_yr
         239             :     case( 'SERIAL' )
         240             :     case default
         241           0 :        write(iulog,*) 'trcdata_init: invalid data type: '//trim(data_type)//' file: '//trim(filename)
         242           0 :        write(iulog,*) 'trcdata_init: valid data types: SERIAL | CYCLICAL | CYCLICAL_LIST | FIXED | INTERP_MISSING_MONTHS '
         243        6144 :        call endrun('trcdata_init: invalid data type: '//trim(data_type)//' file: '//trim(filename))
         244             :     endselect
         245             :
         246        6144 :     if ( (.not.file%fixed) .and. ((data_fixed_ymd>0._r8) .or.(data_fixed_tod>0._r8))) then
         247           0 :        call endrun('trcdata_init: Cannot specify data_fixed_ymd or data_fixed_tod if data type is not FIXED')
         248             :     endif
         249        6144 :     if ( (.not.file%cyclical) .and. (data_cycle_yr>0._r8) ) then
         250           0 :        call endrun('trcdata_init: Cannot specify data_cycle_yr if data type is not CYCLICAL')
         251             :     endif
         252             :
         253        6144 :     if (file%top_bndry .and. file%top_layer) then
         254           0 :        call endrun('trcdata_init: Cannot set both file%top_bndry and file%top_layer to TRUE.')
         255             :     end if
         256             :
         257        6144 :     if (masterproc) then
         258           8 :        write(iulog,*) 'trcdata_init: data type: '//trim(data_type)//' file: '//trim(filename)
         259             :     endif
         260             :
         261             :     ! if there is no list of files (len_trim(file%filenames_list)<1) then
         262             :     !  -> set curr_filename from namelist rather from restart data
         263        6144 :     if ( len_trim(file%curr_filename)<1 .or. len_trim(file%filenames_list)<1 .or. file%fixed ) then ! initial run
         264        6144 :        file%curr_filename = trim(filename)
         265             :
         266        6144 :        call get_model_time(file)
         267             :
         268        6144 :        if ( file%fixed ) then
         269           0 :           str_yr = data_fixed_ymd/10000
         270           0 :           str_mon = (data_fixed_ymd - str_yr*10000)/100
         271           0 :           str_day = data_fixed_ymd - str_yr*10000 - str_mon*100
         272           0 :           call set_time_float_from_date( start_time, str_yr, str_mon, str_day, data_fixed_tod )
         273           0 :           file%offset_time = start_time - file%curr_mod_time
         274             :        else
         275        6144 :           file%offset_time = 0
         276             :        endif
         277             :     endif
         278             :
         279        6144 :     call set_time_float_from_date( time2, 2, 1, 1, 0 )
         280        6144 :     call set_time_float_from_date( time1, 1, 1, 1, 0 )
         281        6144 :     file%one_yr = time2-time1
         282             :
         283        6144 :     if ( file%cyclical .or. file%cyclical_list) then
         284        3072 :        file%cyc_ndx_beg = -1
         285        3072 :        file%cyc_ndx_end = -1
         286        3072 :        if ( file%cyc_yr /= 0 ) then
         287        3072 :           call set_time_float_from_date( time1, file%cyc_yr  , 1, 1, 0 )
         288        3072 :           call set_time_float_from_date( time2, file%cyc_yr+1, 1, 1, 0 )
         289        3072 :           file%one_yr = time2-time1
         290             :        endif
         291             :
         292        3072 :        if ( file%cyclical ) then
         293             :           call open_trc_datafile( file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times, &
         294        3072 :                cyc_ndx_beg=file%cyc_ndx_beg, cyc_ndx_end=file%cyc_ndx_end, cyc_yr=file%cyc_yr )
         295             :        else
         296           0 :           call open_trc_datafile( file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times )
         297             :        endif
         298             :     else
         299        3072 :        call open_trc_datafile( file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times )
         300    37011456 :        file%curr_data_times = file%curr_data_times - file%offset_time
         301             :     endif
         302             :
         303        6144 :     call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling)
         304        6144 :     ierr = pio_inq_dimid( file%curr_fileid, 'ncol', idx )
         305        6144 :     file%unstructured = (ierr==PIO_NOERR)
         306        6144 :     if (.not.file%unstructured) then
         307        6144 :        ierr = pio_inq_dimid( file%curr_fileid, 'lon', idx )
         308        6144 :        file%zonal_ave = (ierr/=PIO_NOERR)
         309             :     endif
         310        6144 :     call pio_seterrorhandling(File%curr_fileid, err_handling)
         311             :
         312        6144 :     plon = get_dyn_grid_parm('plon')
         313        6144 :     plat = get_dyn_grid_parm('plat')
         314             :
         315        6144 :     if ( file%zonal_ave ) then
         316             :
         317        3072 :        file%nlon = 1
         318             :
         319             :     else
         320             :
         321        3072 :        if (.not. file%unstructured ) then
         322        3072 :           call get_dimension( file%curr_fileid, 'lon', file%nlon, dimid=old_dimid, data=file%lons )
         323             :
         324      445440 :           file%lons =  file%lons * d2r
         325             :
         326        3072 :           lon_dimid = old_dimid
         327             :        end if
         328             :     endif
         329             :
         330        6144 :     ierr = pio_inq_dimid( file%curr_fileid, 'time', old_dimid)
         331             :
         332        6144 :     if (.not. file%unstructured ) then
         333             :        ! Hack to work with weird netCDF and old gcc or NAG bug.
         334        6144 :        tim_dimid = old_dimid
         335             :
         336        6144 :        call get_dimension( file%curr_fileid, 'lat', file%nlat, dimid=old_dimid, data=file%lats )
         337      890880 :        file%lats =  file%lats * d2r
         338             :
         339        6144 :        lat_dimid = old_dimid
         340             :     endif
         341             :
         342        6144 :     call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling)
         343        6144 :     ierr = pio_inq_varid( file%curr_fileid, 'PS', file%ps_id )
         344        6144 :     file%has_ps = (ierr==PIO_NOERR)
         345        6144 :     ierr = pio_inq_dimid( file%curr_fileid, 'altitude', idx )
         346        6144 :     file%alt_data = (ierr==PIO_NOERR)
         347             :
         348        6144 :     call pio_seterrorhandling(File%curr_fileid, err_handling)
         349             :
         350        6144 :     if ( file%has_ps .and. .not.file%unstructured ) then
         351        4608 :        if ( file%zonal_ave ) then
         352        3072 :           ierr = pio_inq_vardimid (file%curr_fileid, file%ps_id, dimids(1:2))
         353        9216 :           do did = 1,2
         354        9216 :              if      ( dimids(did) == lat_dimid ) then
         355        3072 :                 file%ps_coords(LATDIM) = did
         356        3072 :                 file%ps_order(did) = LATDIM
         357        3072 :              else if ( dimids(did) == tim_dimid ) then
         358        3072 :                 file%ps_coords(PS_TIMDIM) = did
         359        3072 :                 file%ps_order(did) = PS_TIMDIM
         360             :              endif
         361             :           enddo
         362             :        else
         363        1536 :           ierr = pio_inq_vardimid (file%curr_fileid, file%ps_id, dimids(1:3))
         364        6144 :           do did = 1,3
         365        6144 :              if      ( dimids(did) == lon_dimid ) then
         366        1536 :                 file%ps_coords(LONDIM) = did
         367        1536 :                 file%ps_order(did) = LONDIM
         368        3072 :              else if ( dimids(did) == lat_dimid ) then
         369        1536 :                 file%ps_coords(LATDIM) = did
         370        1536 :                 file%ps_order(did) = LATDIM
         371        1536 :              else if ( dimids(did) == tim_dimid ) then
         372        1536 :                 file%ps_coords(PS_TIMDIM) = did
         373        1536 :                 file%ps_order(did) = PS_TIMDIM
         374             :              endif
         375             :           enddo
         376             :        end if
         377             :     endif
         378             :
         379        6144 :     if (masterproc) then
         380           8 :        write(iulog,*) 'trcdata_init: file%has_ps = ' , file%has_ps
         381             :     endif ! masterproc
         382             :
         383        6144 :     if (file%alt_data) then
         384           0 :        call get_dimension( file%curr_fileid, 'altitude_int', file%nilev,  data=file%ilevs  )
         385           0 :        call get_dimension( file%curr_fileid, 'altitude',     file%nlev, dimid=old_dimid, data=file%levs  )
         386             :     else
         387        6144 :        call get_dimension( file%curr_fileid, 'lev', file%nlev, dimid=old_dimid, data=file%levs  )
         388        6144 :        if (old_dimid>0) then
         389      259584 :           file%levs =  file%levs*100._r8 ! mbar->pascals
         390             :        endif
         391             :     endif
         392             :
         393             :     ! For some bizarre reason, netCDF with older gcc is keeping a pointer to the dimid, and overwriting it later!
         394             :     ! Hackish workaround is to make a copy...
         395        6144 :     lev_dimid = old_dimid
         396             :
         397        6144 :     if (file%has_ps) then
         398             :
         399       18432 :        allocate( file%hyam(file%nlev),  file%hybm(file%nlev), stat=astat )
         400        4608 :        if( astat /= 0 ) then
         401           0 :           write(iulog,*) 'trcdata_init: file%hyam,file%hybm allocation error = ',astat
         402           0 :           call endrun('trcdata_init: failed to allocate file%hyam and file%hybm arrays')
         403             :        end if
         404             :
         405       18432 :        allocate( file%hyai(file%nlev+1),  file%hybi(file%nlev+1), stat=astat )
         406        4608 :        if( astat /= 0 ) then
         407           0 :           write(iulog,*) 'trcdata_init: file%hyai,file%hybi allocation error = ',astat
         408           0 :           call endrun('trcdata_init: failed to allocate file%hyai and file%hybi arrays')
         409             :        end if
         410             :
         411        4608 :        call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling)
         412        4608 :        ierr = pio_inq_varid( file%curr_fileid, 'P0', varid)
         413        4608 :        call pio_seterrorhandling(File%curr_fileid, err_handling)
         414             :
         415        4608 :        if ( ierr == PIO_NOERR ) then
         416        4608 :           ierr = pio_get_var( file%curr_fileid, varid, file%p0 )
         417             :        else
         418           0 :           file%p0 = 100000._r8
         419             :        endif
         420        4608 :        ierr = pio_inq_varid( file%curr_fileid, 'hyam', varid )
         421        4608 :        ierr = pio_get_var( file%curr_fileid, varid, file%hyam )
         422        4608 :        ierr = pio_inq_varid( file%curr_fileid, 'hybm', varid )
         423        4608 :        ierr = pio_get_var( file%curr_fileid, varid, file%hybm )
         424        4608 :        if (file%conserve_column) then
         425           0 :           ierr = pio_inq_varid( file%curr_fileid, 'hyai', varid )
         426           0 :           ierr = pio_get_var( file%curr_fileid, varid, file%hyai )
         427           0 :           ierr = pio_inq_varid( file%curr_fileid, 'hybi', varid )
         428           0 :           ierr = pio_get_var( file%curr_fileid, varid, file%hybi )
         429             :        endif
         430             :
         431       13824 :        allocate( file%ps_in(1)%data(pcols,begchunk:endchunk), stat=astat   )
         432        4608 :        if( astat/= 0 ) then
         433           0 :           write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(1)%data array; error = ',astat
         434           0 :           call endrun
         435             :        end if
         436       13824 :        allocate( file%ps_in(2)%data(pcols,begchunk:endchunk), stat=astat   )
         437        4608 :        if( astat/= 0 ) then
         438           0 :           write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(2)%data array; error = ',astat
         439           0 :           call endrun
         440             :        end if
         441        4608 :        if( file%fill_in_months ) then
         442           0 :           allocate( file%ps_in(3)%data(pcols,begchunk:endchunk), stat=astat   )
         443           0 :           if( astat/= 0 ) then
         444           0 :              write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(3)%data array; error = ',astat
         445           0 :              call endrun
         446             :           end if
         447           0 :           allocate( file%ps_in(4)%data(pcols,begchunk:endchunk), stat=astat   )
         448           0 :           if( astat/= 0 ) then
         449           0 :              write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(4)%data array; error = ',astat
         450           0 :              call endrun
         451             :           end if
         452             :        end if
         453             :     endif
         454             :
         455             :
         456        6144 :     call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling)
         457             :
         458       59904 :     flds_loop: do f = 1,mxnflds
         459             :
         460             :        ! get netcdf variable id for the field
         461       53760 :        ierr = pio_inq_varid( file%curr_fileid, flds(f)%srcnam, flds(f)%var_id )
         462       53760 :        if (ierr/=pio_noerr) then
         463           0 :           call endrun('trcdata_init: Cannot find var "'//trim(flds(f)%srcnam)// &
         464           0 :                       '" in file "'//trim(file%curr_filename)//'"')
         465             :        endif
         466             :
         467             :        ! determine if the field has a vertical dimension
         468             :
         469       53760 :        if (lev_dimid>0) then
         470       32256 :           ierr = pio_inquire_variable(  file%curr_fileid, flds(f)%var_id, ndims=nvardims )
         471       32256 :           ierr = pio_inquire_variable(  file%curr_fileid, flds(f)%var_id, dimids=vardimids(:nvardims) )
         472       84480 :           flds(f)%srf_fld = .not.any(vardimids(:nvardims)==lev_dimid)
         473             :        else
         474       21504 :           flds(f)%srf_fld = .true.
         475             :        endif
         476             :
         477             :        ! allocate memory only if not already in pbuf2d
         478             :
         479       53760 :        if ( .not. file%in_pbuf(f) ) then
         480       21504 :           if ( flds(f)%srf_fld .or. file%top_bndry .or. file%top_layer ) then
         481       64512 :              allocate( flds(f)%data(pcols,1,begchunk:endchunk), stat=astat   )
         482             :           else
         483           0 :              allocate( flds(f)%data(pcols,pver,begchunk:endchunk), stat=astat   )
         484             :           endif
         485       21504 :           if( astat/= 0 ) then
         486           0 :              write(iulog,*) 'trcdata_init: failed to allocate flds(f)%data array; error = ',astat
         487           0 :              call endrun
         488             :           end if
         489             :        else
         490       32256 :           flds(f)%pbuf_ndx = pbuf_get_index(flds(f)%fldnam,errcode)
         491             :        endif
         492             :
         493       53760 :        if (flds(f)%srf_fld) then
         494       64512 :           allocate( flds(f)%input(1)%data(pcols,1,begchunk:endchunk), stat=astat   )
         495             :        else
         496      129024 :           allocate( flds(f)%input(1)%data(pcols,file%nlev,begchunk:endchunk), stat=astat   )
         497             :        endif
         498       53760 :        if( astat/= 0 ) then
         499           0 :           write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(1)%data array; error = ',astat
         500           0 :           call endrun
         501             :        end if
         502       53760 :        if (flds(f)%srf_fld) then
         503       64512 :           allocate( flds(f)%input(2)%data(pcols,1,begchunk:endchunk), stat=astat   )
         504             :        else
         505      129024 :           allocate( flds(f)%input(2)%data(pcols,file%nlev,begchunk:endchunk), stat=astat   )
         506             :        endif
         507       53760 :        if( astat/= 0 ) then
         508           0 :           write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(2)%data array; error = ',astat
         509           0 :           call endrun
         510             :        end if
         511             :
         512       53760 :        if( file%fill_in_months ) then
         513           0 :           if (flds(f)%srf_fld) then
         514           0 :              allocate( flds(f)%input(3)%data(pcols,1,begchunk:endchunk), stat=astat   )
         515             :           else
         516           0 :              allocate( flds(f)%input(3)%data(pcols,file%nlev,begchunk:endchunk), stat=astat   )
         517             :           endif
         518           0 :           if( astat/= 0 ) then
         519           0 :              write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(3)%data array; error = ',astat
         520           0 :              call endrun
         521             :           end if
         522           0 :           if (flds(f)%srf_fld) then
         523           0 :              allocate( flds(f)%input(4)%data(pcols,1,begchunk:endchunk), stat=astat   )
         524             :           else
         525           0 :              allocate( flds(f)%input(4)%data(pcols,file%nlev,begchunk:endchunk), stat=astat   )
         526             :           endif
         527           0 :           if( astat/= 0 ) then
         528           0 :              write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(4)%data array; error = ',astat
         529           0 :              call endrun
         530             :           end if
         531             :        endif
         532             :
         533       53760 :        if ( file%zonal_ave ) then
         534       12288 :           ierr = pio_inq_vardimid (file%curr_fileid, flds(f)%var_id, dimids(1:3))
         535       49152 :           do did = 1,3
         536       49152 :              if      ( dimids(did) == lat_dimid ) then
         537       12288 :                 flds(f)%coords(ZA_LATDIM) = did
         538       12288 :                 flds(f)%order(did) = ZA_LATDIM
         539       24576 :              else if ( dimids(did) == lev_dimid ) then
         540       12288 :                 flds(f)%coords(ZA_LEVDIM) = did
         541       12288 :                 flds(f)%order(did) = ZA_LEVDIM
         542       12288 :              else if ( dimids(did) == tim_dimid ) then
         543       12288 :                 flds(f)%coords(ZA_TIMDIM) = did
         544       12288 :                 flds(f)%order(did) = ZA_TIMDIM
         545             :              endif
         546             :           enddo
         547       41472 :        else if ( flds(f)%srf_fld .and. .not.file%unstructured ) then
         548       21504 :           ierr = pio_inq_vardimid (file%curr_fileid, flds(f)%var_id, dimids(1:3))
         549       86016 :           do did = 1,3
         550       86016 :              if      ( dimids(did) == lon_dimid ) then
         551       21504 :                 flds(f)%coords(LONDIM) = did
         552       21504 :                 flds(f)%order(did) = LONDIM
         553       43008 :              else if ( dimids(did) == lat_dimid ) then
         554       21504 :                 flds(f)%coords(LATDIM) = did
         555       21504 :                 flds(f)%order(did) = LATDIM
         556       21504 :              else if ( dimids(did) == tim_dimid ) then
         557       21504 :                 flds(f)%coords(PS_TIMDIM) = did
         558       21504 :                 flds(f)%order(did) = PS_TIMDIM
         559             :              endif
         560             :           enddo
         561       19968 :        else if (.not. file%unstructured ) then
         562       19968 :           ierr = pio_inq_vardimid (file%curr_fileid, flds(f)%var_id, dimids)
         563       99840 :           do did = 1,4
         564       99840 :              if      ( dimids(did) == lon_dimid ) then
         565       19968 :                 flds(f)%coords(LONDIM) = did
         566       19968 :                 flds(f)%order(did) = LONDIM
         567       59904 :              else if ( dimids(did) == lat_dimid ) then
         568       19968 :                 flds(f)%coords(LATDIM) = did
         569       19968 :                 flds(f)%order(did) = LATDIM
         570       39936 :              else if ( dimids(did) == lev_dimid ) then
         571       19968 :                 flds(f)%coords(LEVDIM) = did
         572       19968 :                 flds(f)%order(did) = LEVDIM
         573       19968 :              else if ( dimids(did) == tim_dimid ) then
         574       19968 :                 flds(f)%coords(TIMDIM) = did
         575       19968 :                 flds(f)%order(did) = TIMDIM
         576             :              endif
         577             :           enddo
         578             :        endif
         579             :
         580       53760 :        ierr = pio_get_att( file%curr_fileid, flds(f)%var_id, 'units', data_units)
         581       59904 :        flds(f)%units = trim(data_units(1:32))
         582             :
         583             :     enddo flds_loop
         584             :
         585        6144 :     call pio_seterrorhandling(File%curr_fileid, err_handling)
         586             :
         587             : ! if weighting by latitude, compute weighting for horizontal interpolation
         588        6144 :     if( file%weight_by_lat ) then
         589           0 :         if(dycore_is('UNSTRUCTURED') ) then
         590           0 :           call endrun('trcdata_init: weighting by latitude not implemented for unstructured grids')
         591             :         endif
         592             :
         593             : ! get dimensions of CAM resolution
         594           0 :         plon = get_dyn_grid_parm('plon')
         595           0 :         plat = get_dyn_grid_parm('plat')
         596             :
         597           0 :         allocate(lam(plon), phi(plat))
         598           0 :         call get_horiz_grid_d(plat, clat_d_out=phi)
         599           0 :         call get_horiz_grid_d(plon, clon_d_out=lam)
         600             :
         601           0 :          if(.not.allocated(lon_global_grid_ndx)) allocate(lon_global_grid_ndx(pcols,begchunk:endchunk))
         602           0 :          if(.not.allocated(lat_global_grid_ndx)) allocate(lat_global_grid_ndx(pcols,begchunk:endchunk))
         603           0 :         lon_global_grid_ndx=-huge(1)
         604           0 :         lat_global_grid_ndx=-huge(1)
         605             :
         606           0 :         do lchnk = begchunk, endchunk
         607           0 :            ncol = get_ncols_p(lchnk)
         608           0 :            call get_rlat_all_p(lchnk, ncol, rlats(:ncol))
         609           0 :            call get_rlon_all_p(lchnk, ncol, rlons(:ncol))
         610           0 :            do icol= 1,ncol
         611           0 :               found=.false.
         612           0 :               find_col: do j = 1,plat
         613           0 :                  do i = 1,plon
         614           0 :                     if (rlats(icol)==phi(j) .and. rlons(icol)==lam(i)) then
         615             :                        found=.true.
         616             :                        exit find_col
         617             :                     endif
         618             :                  enddo
         619             :               enddo find_col
         620             :
         621           0 :               if (.not.found) call endrun('trcdata_init: not able find physics column coordinate')
         622           0 :               lon_global_grid_ndx(icol,lchnk) = i
         623           0 :               lat_global_grid_ndx(icol,lchnk) = j
         624             :            end do
         625             :         enddo
         626             :
         627           0 :         deallocate(phi,lam)
         628             :
         629             : ! weight_x & weight_y are weighting function for x & y interpolation
         630           0 :         allocate(file%weight_x(plon,file%nlon), stat=astat)
         631           0 :         if( astat /= 0 ) then
         632           0 :            write(iulog,*) 'trcdata_init: file%weight_x allocation error = ',astat
         633           0 :            call endrun('trcdata_init: failed to allocate weight_x array')
         634             :         end if
         635           0 :         allocate(file%weight_y(plat,file%nlat), stat=astat)
         636           0 :         if( astat /= 0 ) then
         637           0 :            write(iulog,*) 'trcdata_init: file%weight_y allocation error = ',astat
         638           0 :            call endrun('trcdata_init: failed to allocate weight_y array')
         639             :         end if
         640           0 :         allocate(file%count_x(plon), stat=astat)
         641           0 :         if( astat /= 0 ) then
         642           0 :            write(iulog,*) 'trcdata_init: file%count_x allocation error = ',astat
         643           0 :            call endrun('trcdata_init: failed to allocate count_x array')
         644             :         end if
         645           0 :         allocate(file%count_y(plat), stat=astat)
         646           0 :         if( astat /= 0 ) then
         647           0 :            write(iulog,*) 'trcdata_init: file%count_y allocation error = ',astat
         648           0 :            call endrun('trcdata_init: failed to allocate count_y array')
         649             :         end if
         650           0 :         allocate(file%index_x(plon,file%nlon), stat=astat)
         651           0 :         if( astat /= 0 ) then
         652           0 :            write(iulog,*) 'trcdata_init: file%index_x allocation error = ',astat
         653           0 :            call endrun('trcdata_init: failed to allocate index_x array')
         654             :         end if
         655           0 :         allocate(file%index_y(plat,file%nlat), stat=astat)
         656           0 :         if( astat /= 0 ) then
         657           0 :            write(iulog,*) 'trcdata_init: file%index_y allocation error = ',astat
         658           0 :            call endrun('trcdata_init: failed to allocate index_y array')
         659             :         end if
         660           0 :         file%weight_x(:,:) = 0.0_r8
         661           0 :         file%weight_y(:,:) = 0.0_r8
         662           0 :         file%count_x(:) = 0
         663           0 :         file%count_y(:) = 0
         664           0 :         file%index_x(:,:) = 0
         665           0 :         file%index_y(:,:) = 0
         666             :
         667           0 :         if( file%dist ) then
         668           0 :            allocate(file%weight0_x(plon,file%nlon), stat=astat)
         669           0 :            if( astat /= 0 ) then
         670           0 :               write(iulog,*) 'trcdata_init: file%weight0_x allocation error = ',astat
         671           0 :               call endrun('trcdata_init: failed to allocate weight0_x array')
         672             :            end if
         673           0 :            allocate(file%weight0_y(plat,file%nlat), stat=astat)
         674           0 :            if( astat /= 0 ) then
         675           0 :               write(iulog,*) 'trcdata_init: file%weight0_y allocation error = ',astat
         676           0 :               call endrun('trcdata_init: failed to allocate weight0_y array')
         677             :            end if
         678           0 :            allocate(file%count0_x(plon), stat=astat)
         679           0 :            if( astat /= 0 ) then
         680           0 :               write(iulog,*) 'trcdata_init: file%count0_x allocation error = ',astat
         681           0 :               call endrun('trcdata_init: failed to allocate count0_x array')
         682             :            end if
         683           0 :            allocate(file%count0_y(plat), stat=astat)
         684           0 :            if( astat /= 0 ) then
         685           0 :               write(iulog,*) 'trcdata_init: file%count0_y allocation error = ',astat
         686           0 :               call endrun('trcdata_init: failed to allocate count0_y array')
         687             :            end if
         688           0 :            allocate(file%index0_x(plon,file%nlon), stat=astat)
         689           0 :            if( astat /= 0 ) then
         690           0 :               write(iulog,*) 'trcdata_init: file%index0_x allocation error = ',astat
         691           0 :               call endrun('trcdata_init: failed to allocate index0_x array')
         692             :            end if
         693           0 :            allocate(file%index0_y(plat,file%nlat), stat=astat)
         694           0 :            if( astat /= 0 ) then
         695           0 :               write(iulog,*) 'trcdata_init: file%index0_y allocation error = ',astat
         696           0 :               call endrun('trcdata_init: failed to allocate index0_y array')
         697             :            end if
         698           0 :            file%weight0_x(:,:) = 0.0_r8
         699           0 :            file%weight0_y(:,:) = 0.0_r8
         700           0 :            file%count0_x(:) = 0
         701           0 :            file%count0_y(:) = 0
         702           0 :            file%index0_x(:,:) = 0
         703           0 :            file%index0_y(:,:) = 0
         704             :         endif
         705             :
         706           0 :         if(masterproc) then
         707             : ! compute weighting
         708             :             call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats, &
         709           0 :                                 plon,plat,file%weight_x,file%weight_y,file%dist)
         710             :
         711           0 :             do i2=1,plon
         712           0 :                file%count_x(i2) = 0
         713           0 :                do i1=1,file%nlon
         714           0 :                   if(file%weight_x(i2,i1)>0.0_r8 ) then
         715           0 :                      file%count_x(i2) = file%count_x(i2) + 1
         716           0 :                      file%index_x(i2,file%count_x(i2)) = i1
         717             :                   endif
         718             :                enddo
         719             :             enddo
         720             :
         721           0 :             do j2=1,plat
         722           0 :                file%count_y(j2) = 0
         723           0 :                do j1=1,file%nlat
         724           0 :                   if(file%weight_y(j2,j1)>0.0_r8 ) then
         725           0 :                      file%count_y(j2) = file%count_y(j2) + 1
         726           0 :                      file%index_y(j2,file%count_y(j2)) = j1
         727             :                   endif
         728             :                enddo
         729             :             enddo
         730             :
         731           0 :            if( file%dist ) then
         732             :             call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,&
         733           0 :                                 plon,plat,file%weight0_x,file%weight0_y,file%dist)
         734             :
         735           0 :             do i2=1,plon
         736           0 :                file%count0_x(i2) = 0
         737           0 :                do i1=1,file%nlon
         738           0 :                   if(file%weight0_x(i2,i1)>0.0_r8 ) then
         739           0 :                      file%count0_x(i2) = file%count0_x(i2) + 1
         740           0 :                      file%index0_x(i2,file%count0_x(i2)) = i1
         741             :                   endif
         742             :                enddo
         743             :             enddo
         744             :
         745           0 :             do j2=1,plat
         746           0 :                file%count0_y(j2) = 0
         747           0 :                do j1=1,file%nlat
         748           0 :                   if(file%weight0_y(j2,j1)>0.0_r8 ) then
         749           0 :                      file%count0_y(j2) = file%count0_y(j2) + 1
         750           0 :                      file%index0_y(j2,file%count0_y(j2)) = j1
         751             :                   endif
         752             :                enddo
         753             :             enddo
         754             :            endif
         755             :         endif
         756             :
         757           0 :         call mpi_bcast(file%weight_x, plon*file%nlon, mpi_real8 , mstrid, mpicom,ierr)
         758           0 :         if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_x")
         759           0 :         call mpi_bcast(file%weight_y, plat*file%nlat, mpi_real8 , mstrid, mpicom,ierr)
         760           0 :         if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_y")
         761           0 :         call mpi_bcast(file%count_x, plon, mpi_integer , mstrid, mpicom,ierr)
         762           0 :         if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_x")
         763           0 :         call mpi_bcast(file%count_y, plat, mpi_integer , mstrid, mpicom,ierr)
         764           0 :         if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_y")
         765           0 :         call mpi_bcast(file%index_x, plon*file%nlon, mpi_integer , mstrid, mpicom,ierr)
         766           0 :         if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index_x")
         767           0 :         call mpi_bcast(file%index_y, plat*file%nlat, mpi_integer , mstrid, mpicom,ierr)
         768           0 :         if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index_y")
         769           0 :         if( file%dist ) then
         770           0 :            call mpi_bcast(file%weight0_x, plon*file%nlon, mpi_real8 , mstrid, mpicom,ierr)
         771           0 :            if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_x")
         772           0 :            call mpi_bcast(file%weight0_y,  plat*file%nlat, mpi_real8 , mstrid, mpicom,ierr)
         773           0 :            if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_y")
         774           0 :            call mpi_bcast(file%count0_x, plon, mpi_integer , mstrid, mpicom,ierr)
         775           0 :            if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_x")
         776           0 :            call mpi_bcast(file%count0_y, plon, mpi_integer , mstrid, mpicom,ierr)
         777           0 :            if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_y")
         778           0 :            call mpi_bcast(file%index0_x, plon*file%nlon, mpi_integer , mstrid, mpicom,ierr)
         779           0 :            if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_x")
         780           0 :            call mpi_bcast(file%index0_y,  plat*file%nlat, mpi_integer , mstrid, mpicom,ierr)
         781           0 :            if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_y")
         782             :         endif
         783             :     endif
         784             :
         785        6144 :   end subroutine trcdata_init
         786             :
         787             : !-----------------------------------------------------------------------
         788             : ! Reads more data if needed and interpolates data to current model time
         789             : !-----------------------------------------------------------------------
         790     1483776 :   subroutine advance_trcdata( flds, file, state, pbuf2d )
         791        6144 :     use physics_types,only : physics_state
         792             :
         793             :     implicit none
         794             :
         795             :     type(trfile),        intent(inout) :: file
         796             :     type(trfld),         intent(inout) :: flds(:)
         797             :     type(physics_state), intent(in)    :: state(begchunk:endchunk)
         798             :
         799             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         800             :
         801             :     real(r8) :: data_time
         802             :
         803     1483776 :     call t_startf('advance_trcdata')
         804     1483776 :     if ( .not.( file%fixed .and. file%initialized ) ) then
         805             :
         806     1483776 :        call get_model_time(file)
         807             :
         808     1483776 :        data_time = file%datatimep
         809             :
         810     1483776 :        if ( file%cyclical .or. file%cyclical_list ) then
         811             :           ! wrap around
         812      741888 :           if ( (file%datatimep<file%datatimem) .and. (file%curr_mod_time>file%datatimem) ) then
         813           0 :              data_time = data_time + file%one_yr
         814             :           endif
         815             :        endif
         816             :
         817             :     ! For stepTime need to advance if the times are equal
         818             :     ! Should not impact other runs?
         819     1483776 :        if ( file%curr_mod_time >= data_time ) then
         820        9216 :           call t_startf('read_next_trcdata')
         821        9216 :           call read_next_trcdata( flds, file )
         822        9216 :           call t_stopf('read_next_trcdata')
         823        9302 :           if(masterproc) write(iulog,*) 'READ_NEXT_TRCDATA ', flds%fldnam
         824             :        end if
         825             :
         826             :     endif
         827             :
         828             :     ! need to interpolate the data, regardless
         829             :     ! each mpi task needs to interpolate
         830     1483776 :     call t_startf('interpolate_trcdata')
         831     1483776 :     call interpolate_trcdata( state, flds, file, pbuf2d )
         832     1483776 :     call t_stopf('interpolate_trcdata')
         833             :
         834     1483776 :     file%initialized = .true.
         835             :
         836     1483776 :     call t_stopf('advance_trcdata')
         837             :
         838     1483776 :   end subroutine advance_trcdata
         839             :
         840             : !-------------------------------------------------------------------
         841             : !-------------------------------------------------------------------
         842           0 :   subroutine get_fld_data( flds, field_name, data, ncol, lchnk, pbuf )
         843             :
         844             :
         845             :     implicit none
         846             :
         847             :     type(trfld), intent(inout) :: flds(:)
         848             :     character(len=*), intent(in) :: field_name
         849             :     real(r8), intent(out) :: data(:,:)
         850             :     integer, intent(in) :: lchnk
         851             :     integer, intent(in) :: ncol
         852             :     type(physics_buffer_desc), pointer :: pbuf(:)
         853             :
         854             :
         855             :     integer :: f, nflds
         856           0 :     real(r8),pointer  :: tmpptr(:,:)
         857             :
         858           0 :     data(:,:) = 0._r8
         859           0 :     nflds = size(flds)
         860             :
         861           0 :     do f = 1, nflds
         862           0 :        if ( trim(flds(f)%fldnam) == trim(field_name) ) then
         863           0 :           if ( flds(f)%pbuf_ndx>0 ) then
         864             :              call pbuf_get_field(pbuf, flds(f)%pbuf_ndx, tmpptr)
         865           0 :              data(:ncol,:) = tmpptr(:ncol,:)
         866             :           else
         867           0 :              data(:ncol,:) = flds(f)%data(:ncol,:,lchnk)
         868             :           endif
         869             :        endif
         870             :     enddo
         871             :
         872     1483776 :  end subroutine get_fld_data
         873             :
         874             : !-------------------------------------------------------------------
         875             : !-------------------------------------------------------------------
         876           0 :   subroutine get_fld_ndx( flds, field_name, idx  )
         877             :
         878             :     implicit none
         879             :
         880             :     type(trfld), intent(in) :: flds(:)
         881             :     character(len=*), intent(in) :: field_name
         882             :     integer, intent(out) :: idx
         883             :     integer :: f, nflds
         884             :
         885           0 :     idx = -1
         886           0 :     nflds = size(flds)
         887             :
         888           0 :     do f = 1, nflds
         889           0 :        if ( trim(flds(f)%fldnam) == trim(field_name) ) then
         890           0 :           idx = f
         891           0 :           return
         892             :        endif
         893             :     enddo
         894             :
         895             :   end subroutine get_fld_ndx
         896             :
         897             : !------------------------------------------------------------------------------
         898             : !------------------------------------------------------------------------------
         899     1489920 :   subroutine get_model_time(file)
         900             :     implicit none
         901             :     type(trfile), intent(inout) :: file
         902             :
         903             :     integer yr, mon, day, ncsec  ! components of a date
         904             :
         905     1489920 :     call get_curr_date(yr, mon, day, ncsec)
         906             :
         907     1489920 :     if ( file%cyclical .or. file%cyclical_list) yr = file%cyc_yr
         908     1489920 :     call set_time_float_from_date( file%curr_mod_time, yr, mon, day, ncsec )
         909     1489920 :     file%next_mod_time = file%curr_mod_time + get_step_size()/86400._r8
         910             :
         911     1489920 :   end subroutine get_model_time
         912             :
         913             : !------------------------------------------------------------------------------
         914             : !------------------------------------------------------------------------------
         915           0 :   subroutine check_files( file, fids, itms, times_found)
         916             :
         917             :     implicit none
         918             :
         919             :     type(trfile),      intent(inout) :: file
         920             :     type(file_desc_t), intent(out)   :: fids(2) ! ids of files that contains these recs
         921             :     integer, optional, intent(out)   :: itms(2)
         922             :     logical, optional, intent(inout) :: times_found
         923             :
         924             :     !-----------------------------------------------------------------------
         925             :     !   ... local variables
         926             :     !-----------------------------------------------------------------------
         927             :     logical            :: list_cycled
         928             :
         929           0 :     list_cycled = .false.
         930             :
         931             :    !-----------------------------------------------------------------------
         932             :    !        If next time beyond the end of the time list,
         933             :    !        then increment the filename and move on to the next file
         934             :    !-----------------------------------------------------------------------
         935           0 :     if ((file%next_mod_time > file%curr_data_times(size(file%curr_data_times))).or.file%cyclical_list) then
         936           0 :        if (file%cyclical_list) then
         937           0 :           if ( associated(file%next_data_times) ) then
         938           0 :               if ((file%curr_mod_time > file%datatimep)) then
         939             :
         940           0 :                call advance_file(file)
         941             :
         942             :             endif
         943             :          endif
         944             :
         945             :        endif
         946           0 :        if ( .not. associated(file%next_data_times) ) then
         947             :           ! open next file if not already opened...
         948           0 :           if (file%cyclical_list) then
         949             :               file%next_filename = incr_filename( file%curr_filename, filenames_list=file%filenames_list, datapath=file%pathname ,&
         950           0 :                                                   cyclical_list=file%cyclical_list, list_cycled=list_cycled)
         951             :           else
         952           0 :               file%next_filename = incr_filename( file%curr_filename, filenames_list=file%filenames_list, datapath=file%pathname)
         953             :           endif
         954           0 :           call open_trc_datafile( file%next_filename, file%pathname, file%next_fileid, file%next_data_times )
         955           0 :           file%next_data_times = file%next_data_times - file%offset_time
         956             :        endif
         957             :     endif
         958             :
         959             :     !-----------------------------------------------------------------------
         960             :     !        If using next_data_times and the current is greater than or equal to the next, then
         961             :     !        close the current file, and set up for next file.
         962             :     !-----------------------------------------------------------------------
         963           0 :     if ( associated(file%next_data_times) ) then
         964           0 :        if (file%cyclical_list .and. list_cycled) then    ! special case - list cycled
         965             :
         966           0 :           file%datatimem = file%curr_data_times(size(file%curr_data_times))
         967           0 :           itms(1)=size(file%curr_data_times)
         968           0 :           fids(1)=file%curr_fileid
         969             :
         970           0 :           file%datatimep = file%next_data_times(1)
         971           0 :           itms(2)=1
         972           0 :           fids(2) = file%next_fileid
         973             :
         974           0 :           times_found = .true.
         975             :
         976           0 :        else if (file%curr_mod_time >= file%next_data_times(1)) then
         977             :
         978           0 :           call advance_file(file)
         979             :
         980             :        endif
         981             :     endif
         982             :
         983           0 :   end subroutine check_files
         984             :
         985             : !-----------------------------------------------------------------------
         986             : !-----------------------------------------------------------------------
         987           0 :   function incr_filename( filename, filenames_list, datapath, cyclical_list, list_cycled, abort )
         988             :
         989             :     !-----------------------------------------------------------------------
         990             :     !   ... Increment or decrement a date string withing a filename
         991             :     !           the filename date section is assumed to be of the form
         992             :     !           yyyy-dd-mm
         993             :     !-----------------------------------------------------------------------
         994             :
         995             :     use string_utils, only : incstr
         996             :     use shr_file_mod, only : shr_file_getunit, shr_file_freeunit
         997             :
         998             :     implicit none
         999             :
        1000             :     character(len=*),           intent(in)    :: filename ! present dynamical dataset filename
        1001             :     character(len=*), optional, intent(in)    :: filenames_list
        1002             :     character(len=*), optional, intent(in)    :: datapath
        1003             :     logical         , optional, intent(in)    :: cyclical_list  ! If true, allow list to cycle
        1004             :     logical         , optional, intent(out)   :: list_cycled
        1005             :     logical         , optional, intent(in)    :: abort
        1006             :
        1007             :     character(len=shr_kind_cl)                :: incr_filename         ! next filename in the sequence
        1008             :
        1009             :     ! set new next_filename ...
        1010             :
        1011             :     !-----------------------------------------------------------------------
        1012             :     !   ... local variables
        1013             :     !-----------------------------------------------------------------------
        1014             :     integer :: pos, istat
        1015             :     character(len=shr_kind_cl) :: fn_new, line, filepath
        1016             :     integer :: ios,unitnumber
        1017             :     logical :: abort_run
        1018             :
        1019           0 :     if (present(abort)) then
        1020           0 :        abort_run = abort
        1021             :     else
        1022             :        abort_run = .true.
        1023             :     endif
        1024             :
        1025           0 :     if (present(list_cycled)) list_cycled = .false.
        1026             :
        1027           0 :     if (( .not. present(filenames_list)) .or.(len_trim(filenames_list) == 0)) then
        1028             :        !-----------------------------------------------------------------------
        1029             :        !        ... ccm type filename
        1030             :        !-----------------------------------------------------------------------
        1031           0 :        pos = len_trim( filename )
        1032           0 :        fn_new = filename(:pos)
        1033           0 :        if ( masterproc ) write(iulog,*) 'incr_flnm: old filename = ',trim(fn_new)
        1034           0 :        if( fn_new(pos-2:) == '.nc' ) then
        1035           0 :           pos = pos - 3
        1036             :        end if
        1037           0 :        istat = incstr( fn_new(:pos), 1 )
        1038           0 :        if( istat /= 0 ) then
        1039           0 :           write(iulog,*) 'incr_flnm: incstr returned ', istat
        1040           0 :           write(iulog,*) '           while trying to decrement ',trim( fn_new )
        1041           0 :           call endrun
        1042             :        end if
        1043             :
        1044             :     else
        1045             :
        1046             :        !-------------------------------------------------------------------
        1047             :        !  ... open filenames_list
        1048             :        !-------------------------------------------------------------------
        1049           0 :        if ( masterproc ) write(iulog,*) 'incr_flnm: old filename = ',trim(filename)
        1050           0 :        if ( masterproc ) write(iulog,*) 'incr_flnm: open filenames_list : ',trim(filenames_list)
        1051           0 :        unitnumber = shr_file_getUnit()
        1052           0 :        if ( present(datapath) ) then
        1053           0 :          filepath = trim(datapath) //'/'// trim(filenames_list)
        1054             :        else
        1055           0 :          filepath = trim(filenames_list)
        1056             :        endif
        1057             :
        1058           0 :        open( unit=unitnumber, file=filepath, iostat=ios, status="OLD")
        1059           0 :        if (ios /= 0) then
        1060           0 :           call endrun('not able to open file: '//trim(filepath))
        1061             :        endif
        1062             :
        1063             :        !-------------------------------------------------------------------
        1064             :        !  ...  read file names
        1065             :        !-------------------------------------------------------------------
        1066           0 :        read( unit=unitnumber, fmt='(A)', iostat=ios ) line
        1067           0 :        if (ios /= 0) then
        1068           0 :           if (abort_run) then
        1069           0 :              call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
        1070             :           else
        1071           0 :              fn_new = 'NOT_FOUND'
        1072           0 :              incr_filename = trim(fn_new)
        1073           0 :              return
        1074             :           endif
        1075             :        endif
        1076             :
        1077             :        !-------------------------------------------------------------------
        1078             :        !      If current filename is '', then initialize with the first filename read in
        1079             :        !      and skip this section.
        1080             :        !-------------------------------------------------------------------
        1081           0 :        if (filename /= '') then
        1082             :
        1083             :           !-------------------------------------------------------------------
        1084             :           !       otherwise read until find current filename
        1085             :           !-------------------------------------------------------------------
        1086           0 :           do while( trim(line) /= trim(filename) )
        1087           0 :              read( unit=unitnumber, fmt='(A)', iostat=ios ) line
        1088           0 :              if (ios /= 0) then
        1089           0 :                 if (abort_run) then
        1090           0 :                    call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
        1091             :                 else
        1092           0 :                    fn_new = 'NOT_FOUND'
        1093           0 :                    incr_filename = trim(fn_new)
        1094           0 :                    return
        1095             :                 endif
        1096             :              endif
        1097             :           enddo
        1098             :
        1099             :           !-------------------------------------------------------------------
        1100             :           !      Read next filename
        1101             :           !-------------------------------------------------------------------
        1102           0 :           read( unit=unitnumber, fmt='(A)', iostat=ios ) line
        1103             :
        1104             :           !---------------------------------------------------------------------------------
        1105             :           !       If cyclical_list, then an end of file is not an error, but rather
        1106             :           !       a signal to rewind and start over
        1107             :           !---------------------------------------------------------------------------------
        1108             :
        1109           0 :           if (ios /= 0) then
        1110           0 :              if (present(cyclical_list)) then
        1111           0 :                 if (cyclical_list) then
        1112           0 :                    list_cycled=.true.
        1113           0 :                    rewind(unitnumber)
        1114           0 :                    read( unit=unitnumber, fmt='(A)', iostat=ios ) line
        1115             :                      ! Error here should never happen, but check just in case
        1116           0 :                    if (ios /= 0) then
        1117           0 :                       call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
        1118             :                    endif
        1119             :                 else
        1120           0 :                    call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
        1121             :                 endif
        1122             :              else
        1123           0 :                 if (abort_run) then
        1124           0 :                    call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
        1125             :                 else
        1126           0 :                    fn_new = 'NOT_FOUND'
        1127           0 :                    incr_filename = trim(fn_new)
        1128           0 :                    return
        1129             :                 endif
        1130             :              endif
        1131             :           endif
        1132             :
        1133             :        endif
        1134             :
        1135             :        !---------------------------------------------------------------------------------
        1136             :        !     Assign the current filename and close the filelist
        1137             :        !---------------------------------------------------------------------------------
        1138           0 :        fn_new = trim(line)
        1139             :
        1140           0 :        close(unit=unitnumber)
        1141           0 :        call shr_file_freeUnit(unitnumber)
        1142             :     endif
        1143             :
        1144             :     !---------------------------------------------------------------------------------
        1145             :     !      return the current filename
        1146             :     !---------------------------------------------------------------------------------
        1147           0 :     incr_filename = trim(fn_new)
        1148           0 :     if ( masterproc ) write(iulog,*) 'incr_flnm: new filename = ',trim(incr_filename)
        1149             :
        1150             :   end function incr_filename
        1151             :
        1152             : !------------------------------------------------------------------------------
        1153             : !------------------------------------------------------------------------------
        1154       27648 :   subroutine find_times( itms, fids, time, file, datatimem, datatimep, times_found )
        1155             :
        1156             :     use intp_util, only: findplb
        1157             :
        1158             :     implicit none
        1159             :
        1160             :     type(trfile), intent(in) :: file
        1161             :     real(r8), intent(out) :: datatimem, datatimep
        1162             :
        1163             :     integer, intent(out) :: itms(2) ! record numbers that bracket time
        1164             :     type(file_desc_t), intent(out) :: fids(2) ! ids of files that contains these recs
        1165             :
        1166             :     real(r8), intent(in) :: time    ! time of interest
        1167             :     logical, intent(inout)  :: times_found
        1168             :
        1169             :     integer :: np1        ! current forward time index of dataset
        1170             :     integer :: n,i      !
        1171             :     integer :: curr_tsize, next_tsize, all_tsize
        1172             :     integer :: astat
        1173             :     integer :: cyc_tsize
        1174             :
        1175        9216 :     real(r8), allocatable, dimension(:):: all_data_times
        1176             :
        1177        9216 :     curr_tsize = size(file%curr_data_times)
        1178        9216 :     next_tsize = 0
        1179        9216 :     if ( associated(file%next_data_times)) next_tsize = size(file%next_data_times)
        1180             :
        1181        9216 :     all_tsize = curr_tsize + next_tsize
        1182             :
        1183       27648 :     allocate( all_data_times( all_tsize ), stat=astat )
        1184        9216 :     if( astat/= 0 ) then
        1185           0 :        write(iulog,*) 'find_times: failed to allocate all_data_times array; error = ',astat
        1186           0 :        call endrun
        1187             :     end if
        1188             :
        1189    76956672 :     all_data_times(:curr_tsize) = file%curr_data_times(:)
        1190        9216 :     if (next_tsize > 0) all_data_times(curr_tsize+1:all_tsize) = file%next_data_times(:)
        1191             :
        1192        9216 :     if ( .not. file%cyclical ) then
        1193        6144 :        if ( all( all_data_times(:) > time ) ) then
        1194           0 :           write(iulog,*) 'FIND_TIMES: ALL data times are after ', time
        1195           0 :           write(iulog,*) 'FIND_TIMES: file: ', trim(file%curr_filename)
        1196           0 :           write(iulog,*) 'FIND_TIMES: time: ', time
        1197           0 :           call endrun('find_times: all(all_data_times(:) > time) '// trim(file%curr_filename) )
        1198             :        endif
        1199             :
        1200             :        ! find bracketing times
        1201    57870336 :        find_times_loop : do n=1, all_tsize-1
        1202    57870336 :           np1 = n + 1
        1203    57870336 :           datatimem = all_data_times(n)   !+ file%offset_time
        1204    57870336 :           datatimep = all_data_times(np1) !+ file%offset_time
        1205             :        ! When stepTime, datatimep may not equal the time (as only datatimem is used)
        1206             :        ! Should not break other runs?
        1207    57870336 :           if ( (time >= datatimem) .and. (time < datatimep) ) then
        1208        6144 :              times_found = .true.
        1209        6144 :              exit find_times_loop
        1210             :           endif
        1211             :        enddo find_times_loop
        1212             :
        1213             :     else  ! file%cyclical
        1214             :
        1215        3072 :        cyc_tsize = file%cyc_ndx_end - file%cyc_ndx_beg + 1
        1216             :
        1217        3072 :        if ( cyc_tsize > 1 ) then
        1218             :
        1219        3072 :           call findplb(all_data_times(file%cyc_ndx_beg:file%cyc_ndx_end),cyc_tsize, time, n )
        1220             :
        1221        3072 :           if (n == cyc_tsize) then
        1222             :              np1 = 1
        1223             :           else
        1224           0 :              np1 = n+1
        1225             :           endif
        1226             :
        1227        3072 :           datatimem = all_data_times(n  +file%cyc_ndx_beg-1)
        1228        3072 :           datatimep = all_data_times(np1+file%cyc_ndx_beg-1)
        1229        3072 :           times_found = .true.
        1230             :
        1231             :        endif
        1232             :     endif
        1233             :
        1234        9216 :     if ( .not. times_found ) then
        1235           0 :        if (masterproc) then
        1236           0 :           write(iulog,*)'FIND_TIMES: Failed to find dates bracketing desired time =', time
        1237           0 :           write(iulog,*) 'filename = '//trim(file%curr_filename)
        1238           0 :           write(iulog,*)' datatimem = ',file%datatimem
        1239           0 :           write(iulog,*)' datatimep = ',file%datatimep
        1240             :        endif
        1241           0 :        return
        1242             :     endif
        1243             :
        1244        9216 :     deallocate( all_data_times, stat=astat )
        1245        9216 :     if( astat/= 0 ) then
        1246           0 :        write(iulog,*) 'find_times: failed to deallocate all_data_times array; error = ',astat
        1247           0 :        call endrun
        1248             :     end if
        1249             :
        1250        9216 :     if ( .not. file%cyclical ) then
        1251        6144 :       itms(1) = n
        1252        6144 :       itms(2) = np1
        1253             :     else
        1254        3072 :       itms(1) = n   +file%cyc_ndx_beg-1
        1255        3072 :       itms(2) = np1 +file%cyc_ndx_beg-1
        1256             :     endif
        1257             :
        1258       27648 :     fids(:) = file%curr_fileid
        1259             :
        1260       27648 :     do i=1,2
        1261       27648 :        if ( itms(i) > curr_tsize ) then
        1262           0 :           itms(i) = itms(i) - curr_tsize
        1263           0 :           fids(i) = file%next_fileid
        1264             :        endif
        1265             :     enddo
        1266             :
        1267        9216 :   end subroutine find_times
        1268             :
        1269             : !------------------------------------------------------------------------
        1270             : !------------------------------------------------------------------------
        1271        9216 :   subroutine read_next_trcdata( flds, file )
        1272             :     implicit none
        1273             :
        1274             :     type (trfile), intent(inout) :: file
        1275             :     type (trfld),intent(inout) :: flds(:)
        1276             :
        1277             :     integer :: recnos(4),i,f,nflds      !
        1278             :     integer :: cnt4(4)            ! array of counts for each dimension
        1279             :     integer :: strt4(4)           ! array of starting indices
        1280             :     integer :: cnt3(3)            ! array of counts for each dimension
        1281             :     integer :: strt3(3)           ! array of starting indices
        1282       46080 :     type(file_desc_t) :: fids(4)
        1283             :     logical :: times_found
        1284             :
        1285             :     integer :: cur_yr, cur_mon, cur_day, cur_sec, yr1, yr2, mon, date, sec
        1286             :     real(r8) :: series1_time, series2_time
        1287             :     type(file_desc_t) :: fid1, fid2
        1288             :
        1289        9216 :     nflds = size(flds)
        1290        9216 :     times_found = .false.
        1291             :
        1292       18432 :     do while( .not. times_found )
        1293        9216 :        call find_times( recnos, fids, file%curr_mod_time, file,file%datatimem, file%datatimep, times_found )
        1294       18432 :        if ( .not. times_found ) then
        1295           0 :           call check_files( file, fids, recnos, times_found )
        1296             :        endif
        1297             :     enddo
        1298             :
        1299             :     !--------------------------------------------------------------
        1300             :     !       If stepTime, then no time interpolation is to be done
        1301             :     !--------------------------------------------------------------
        1302        9216 :     if (file%stepTime) then
        1303           0 :        file%interp_recs = 1
        1304             :     else
        1305        9216 :        file%interp_recs = 2
        1306             :     end if
        1307             :
        1308        9216 :     if ( file%fill_in_months ) then
        1309             :
        1310           0 :        if( file%datatimep-file%datatimem > file%one_yr ) then
        1311             :
        1312           0 :           call get_curr_date(cur_yr, cur_mon, cur_day, cur_sec)
        1313             :
        1314           0 :           call set_date_from_time_float(file%datatimem, yr1, mon, date, sec )
        1315           0 :           call set_date_from_time_float(file%datatimep, yr2, mon, date, sec )
        1316             :
        1317           0 :           call set_time_float_from_date( series1_time, yr1, cur_mon, cur_day, cur_sec )
        1318           0 :           call set_time_float_from_date( series2_time, yr2, cur_mon, cur_day, cur_sec )
        1319             :
        1320           0 :           fid1 = fids(1)
        1321           0 :           fid2 = fids(2)
        1322           0 :           file%cyclical = .true.
        1323           0 :           call set_cycle_indices( fid1, file%cyc_ndx_beg, file%cyc_ndx_end, yr1)
        1324           0 :           call find_times( recnos(1:2), fids(1:2), series1_time, file, file%datatimes(1), file%datatimes(2), times_found )
        1325             :
        1326           0 :           if ( .not. times_found ) then
        1327           0 :               call endrun('read_next_trcdata: time not found for series1_time')
        1328             :           endif
        1329           0 :           call set_cycle_indices( fid2, file%cyc_ndx_beg, file%cyc_ndx_end, yr2)
        1330             :
        1331           0 :           if ( fid1%fh /= fid2%fh ) then
        1332           0 :             file%cyc_ndx_beg = file%cyc_ndx_beg + size(file%curr_data_times)
        1333           0 :             file%cyc_ndx_end = file%cyc_ndx_end + size(file%curr_data_times)
        1334             :           endif
        1335           0 :           call find_times( recnos(3:4), fids(3:4), series2_time, file, file%datatimes(3), file%datatimes(4), times_found )
        1336           0 :           if ( .not. times_found ) then
        1337           0 :               call endrun('read_next_trcdata: time not found for series2_time')
        1338             :           endif
        1339           0 :           file%cyclical = .false.
        1340           0 :           file%interp_recs = 4
        1341             :
        1342           0 :           call set_date_from_time_float( file%datatimes(1), yr1, mon, date, sec )
        1343           0 :           call set_time_float_from_date( file%datatimem, cur_yr,  mon, date, sec )
        1344           0 :           if (file%datatimes(1) > file%datatimes(2) ) then ! wrap around
        1345           0 :             if ( cur_mon == 1 ) then
        1346           0 :                call set_time_float_from_date( file%datatimem, cur_yr-1,  mon, date, sec )
        1347             :             endif
        1348             :           endif
        1349             :
        1350           0 :           call set_date_from_time_float( file%datatimes(2), yr1, mon, date, sec )
        1351           0 :           call set_time_float_from_date( file%datatimep, cur_yr,  mon, date, sec )
        1352           0 :           if (file%datatimes(1) > file%datatimes(2) ) then ! wrap around
        1353           0 :             if ( cur_mon == 12 ) then
        1354           0 :               call set_time_float_from_date( file%datatimep, cur_yr+1,  mon, date, sec )
        1355             :             endif
        1356             :           endif
        1357             :
        1358             :        endif
        1359             :
        1360             :     endif
        1361             :
        1362             :     !
        1363             :     ! Set up hyperslab corners
        1364             :     !
        1365             :
        1366       27648 :     do i=1,file%interp_recs
        1367             :
        1368       92160 :        strt4(:) = 1
        1369       73728 :        strt3(:) = 1
        1370             :
        1371      150528 :        do f = 1,nflds
        1372      150528 :           if ( file%zonal_ave ) then
        1373       49152 :              cnt3(flds(f)%coords(ZA_LATDIM)) = file%nlat
        1374       49152 :              if (flds(f)%srf_fld) then
        1375           0 :                 cnt3(flds(f)%coords(ZA_LEVDIM)) = 1
        1376             :              else
        1377       49152 :                 cnt3(flds(f)%coords(ZA_LEVDIM)) = file%nlev
        1378             :              endif
        1379       49152 :              cnt3(flds(f)%coords(ZA_TIMDIM)) = 1
        1380       49152 :              strt3(flds(f)%coords(ZA_TIMDIM)) = recnos(i)
        1381             :              call read_za_trc( fids(i), flds(f)%var_id, flds(f)%input(i)%data, strt3, cnt3, file, &
        1382      147456 :                   (/ flds(f)%order(ZA_LATDIM),flds(f)%order(ZA_LEVDIM) /) )
        1383       82944 :           else if ( flds(f)%srf_fld ) then
        1384       43008 :              if ( file%unstructured ) then
        1385             :                 ! read data directly onto the unstructureed phys grid -- assumes input data is on same grid as phys
        1386           0 :                 call read_physgrid_2d( fids(i), flds(f)%fldnam,  recnos(i), flds(f)%input(i)%data(:,1,:) )
        1387             :              else
        1388       43008 :                 cnt3( flds(f)%coords(LONDIM)) = file%nlon
        1389       43008 :                 cnt3( flds(f)%coords(LATDIM)) = file%nlat
        1390       43008 :                 cnt3( flds(f)%coords(PS_TIMDIM)) = 1
        1391       43008 :                 strt3(flds(f)%coords(PS_TIMDIM)) = recnos(i)
        1392             :                 call read_2d_trc( fids(i), flds(f)%var_id, flds(f)%input(i)%data(:,1,:), strt3, cnt3, file, &
        1393      129024 :                      (/ flds(f)%order(LONDIM),flds(f)%order(LATDIM) /) )
        1394             :              endif
        1395             :           else
        1396       39936 :              if ( file%unstructured ) then
        1397             :                 ! read data directly onto the unstructureed phys grid -- assumes input data is on same grid as phys
        1398           0 :                 if ( file%alt_data ) then
        1399           0 :                    call read_physgrid_3d( fids(i), flds(f)%fldnam, 'altitude', file%nlev, recnos(i), flds(f)%input(i)%data(:,:,:) )
        1400             :                 else
        1401           0 :                    call read_physgrid_3d( fids(i), flds(f)%fldnam, 'lev', file%nlev, recnos(i), flds(f)%input(i)%data(:,:,:) )
        1402             :                 end if
        1403             :              else
        1404       39936 :                 cnt4(flds(f)%coords(LONDIM)) = file%nlon
        1405       39936 :                 cnt4(flds(f)%coords(LATDIM)) = file%nlat
        1406       39936 :                 cnt4(flds(f)%coords(LEVDIM)) = file%nlev
        1407       39936 :                 cnt4(flds(f)%coords(TIMDIM)) = 1
        1408       39936 :                 strt4(flds(f)%coords(TIMDIM)) = recnos(i)
        1409             :                 call read_3d_trc( fids(i), flds(f)%var_id, flds(f)%input(i)%data, strt4, cnt4, file, &
        1410      159744 :                      (/ flds(f)%order(LONDIM),flds(f)%order(LATDIM),flds(f)%order(LEVDIM) /))
        1411             :              endif
        1412             :
        1413             :           endif
        1414             :
        1415             :        enddo
        1416             :
        1417       27648 :        if ( file%has_ps ) then
        1418       15360 :           if ( file%unstructured ) then
        1419           0 :              call read_physgrid_2d( fids(i), 'PS',  recnos(i), file%ps_in(i)%data )
        1420             :           else
        1421       61440 :              cnt3 = 1
        1422       61440 :              strt3 = 1
        1423       15360 :              if (.not. file%zonal_ave) then
        1424        3072 :                 cnt3(file%ps_coords(LONDIM)) = file%nlon
        1425             :              end if
        1426       15360 :              cnt3(file%ps_coords(LATDIM)) = file%nlat
        1427       15360 :              cnt3(file%ps_coords(PS_TIMDIM)) = 1
        1428       15360 :              strt3(file%ps_coords(PS_TIMDIM)) = recnos(i)
        1429       15360 :              if (file%zonal_ave) then
        1430             :                 call read_2d_trc( fids(i), file%ps_id, file%ps_in(i)%data, strt3(1:2), cnt3(1:2), file, &
        1431       12288 :                      (/ 1, 2 /) )
        1432             :              else
        1433             :                 call read_2d_trc( fids(i), file%ps_id, file%ps_in(i)%data, strt3, cnt3, file, &
        1434        9216 :                      (/ file%ps_order(LONDIM),file%ps_order(LATDIM) /) )
        1435             :              end if
        1436             :           end if
        1437             :        endif
        1438             :
        1439             :     enddo
        1440             :
        1441        9216 :   end subroutine read_next_trcdata
        1442             :
        1443             : !------------------------------------------------------------------------
        1444             :
        1445       58368 :   subroutine read_2d_trc( fid, vid, loc_arr, strt, cnt, file, order )
        1446             :     use interpolate_data,  only : lininterp_init, lininterp, interp_type, lininterp_finish
        1447             :
        1448             :     use ppgrid,       only: pcols, begchunk, endchunk
        1449             :     use phys_grid,    only: get_ncols_p, get_rlat_all_p, get_rlon_all_p
        1450             :     use mo_constants, only: pi
        1451             :     use dycore,       only: dycore_is
        1452             :     use polar_avg,    only: polar_average
        1453             :     use horizontal_interpolate, only : xy_interp
        1454             :
        1455             :     implicit none
        1456             :     type(file_desc_t), intent(in) :: fid
        1457             :     type(var_desc_t), intent(in) :: vid
        1458             :     integer, intent(in) :: strt(:), cnt(:), order(2)
        1459             :     real(r8),intent(out)  :: loc_arr(:,:)
        1460             :     type (trfile), intent(in) :: file
        1461             :
        1462             :     real(r8) :: to_lats(pcols), to_lons(pcols)
        1463       58368 :     real(r8), allocatable, target :: wrk2d(:,:)
        1464       58368 :     real(r8), pointer :: wrk2d_in(:,:)
        1465             :
        1466             :     integer :: c, ierr, ncols
        1467             :     real(r8), parameter :: zero=0_r8, twopi=2_r8*pi
        1468             :     type(interp_type) :: lon_wgts, lat_wgts
        1469             :     integer :: lons(pcols), lats(pcols)
        1470      116736 :     real(r8) :: file_lats(file%nlat)
        1471             :
        1472       58368 :      nullify(wrk2d_in)
        1473      233472 :      allocate( wrk2d(cnt(1),cnt(2)), stat=ierr )
        1474       58368 :      if( ierr /= 0 ) then
        1475           0 :         write(iulog,*) 'read_2d_trc: wrk2d allocation error = ',ierr
        1476           0 :         call endrun
        1477             :      end if
        1478             :
        1479       58368 :      if(order(1)/=1 .or. order(2)/=2 .or. cnt(1)/=file%nlon .or. cnt(2)/=file%nlat) then
        1480       49152 :         allocate( wrk2d_in(file%nlon, file%nlat), stat=ierr )
        1481       12288 :         if( ierr /= 0 ) then
        1482           0 :            write(iulog,*) 'read_2d_trc: wrk2d_in allocation error = ',ierr
        1483           0 :            call endrun
        1484             :         end if
        1485             :      end if
        1486             :
        1487             :
        1488       58368 :     ierr = pio_get_var( fid, vid, strt, cnt, wrk2d )
        1489       58368 :     if(associated(wrk2d_in)) then
        1490     4780032 :        wrk2d_in = reshape( wrk2d(:,:),(/file%nlon,file%nlat/), order=order )
        1491       12288 :        deallocate(wrk2d)
        1492             :     else
        1493       46080 :        wrk2d_in => wrk2d
        1494             :     end if
        1495             :
        1496             :     ! PGI 13.9 bug workaround.
        1497     6841344 :     file_lats = file%lats
        1498             :
        1499             :     ! For zonal average, only interpolate along latitude.
        1500       58368 :     if (file%zonal_ave) then
        1501             :
        1502       61824 :        do c=begchunk,endchunk
        1503       49536 :           ncols = get_ncols_p(c)
        1504       49536 :           call get_rlat_all_p(c, pcols, to_lats)
        1505             :
        1506       49536 :           call lininterp_init(file_lats, file%nlat, to_lats, ncols, 1, lat_wgts)
        1507             :
        1508       49536 :           call lininterp(wrk2d_in(1,:), file%nlat, loc_arr(1:ncols,c-begchunk+1), ncols, lat_wgts)
        1509             :
        1510       61824 :           call lininterp_finish(lat_wgts)
        1511             :        end do
        1512             :
        1513             :     else
        1514             :        ! if weighting by latitude, the perform horizontal interpolation by using weight_x, weight_y
        1515             :
        1516       46080 :        if(file%weight_by_lat) then
        1517             :
        1518           0 :           call t_startf('xy_interp')
        1519             :
        1520           0 :           do c = begchunk,endchunk
        1521           0 :              ncols = get_ncols_p(c)
        1522           0 :              lons(:ncols) = lon_global_grid_ndx(:ncols,c)
        1523           0 :              lats(:ncols) = lat_global_grid_ndx(:ncols,c)
        1524             :
        1525             :              call xy_interp(file%nlon,file%nlat,1,plon,plat,pcols,ncols, &
        1526           0 :                   file%weight_x,file%weight_y,wrk2d_in,loc_arr(:,c-begchunk+1), &
        1527           0 :                   lons,lats,file%count_x,file%count_y,file%index_x,file%index_y)
        1528             :           enddo
        1529             :
        1530           0 :           call t_stopf('xy_interp')
        1531             :
        1532             :        else
        1533      231840 :           do c=begchunk,endchunk
        1534      185760 :              ncols = get_ncols_p(c)
        1535      185760 :              call get_rlat_all_p(c, pcols, to_lats)
        1536      185760 :              call get_rlon_all_p(c, pcols, to_lons)
        1537             :
        1538      185760 :              call lininterp_init(file%lons, file%nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
        1539      185760 :              call lininterp_init(file%lats, file%nlat, to_lats, ncols, 1, lat_wgts)
        1540             :
        1541      185760 :              call lininterp(wrk2d_in, file%nlon, file%nlat, loc_arr(1:ncols,c-begchunk+1), ncols, lon_wgts, lat_wgts)
        1542             :
        1543      185760 :              call lininterp_finish(lon_wgts)
        1544      231840 :              call lininterp_finish(lat_wgts)
        1545             :           end do
        1546             :        endif
        1547             :
        1548             :     end if
        1549             :
        1550       58368 :     if(allocated(wrk2d)) then
        1551       46080 :        deallocate(wrk2d)
        1552             :     else
        1553       12288 :        deallocate(wrk2d_in)
        1554             :     end if
        1555       58368 :     if(dycore_is('LR')) call polar_average(loc_arr)
        1556      116736 :   end subroutine read_2d_trc
        1557             :
        1558             : !------------------------------------------------------------------------
        1559             :
        1560       49152 :   subroutine read_za_trc( fid, vid, loc_arr, strt, cnt, file, order )
        1561       58368 :     use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
        1562             :     use ppgrid,           only : pcols, begchunk, endchunk
        1563             :     use phys_grid,        only : get_ncols_p, get_rlat_all_p
        1564             :
        1565             :     implicit none
        1566             :     type(file_desc_t), intent(in) :: fid
        1567             :     type(var_desc_t),  intent(in) :: vid
        1568             :     integer,           intent(in) :: strt(:), cnt(:)
        1569             :     integer,           intent(in) :: order(2)
        1570             :     real(r8),          intent(out):: loc_arr(:,:,:)
        1571             :     type (trfile),     intent(in) :: file
        1572             :
        1573             :     type(interp_type) :: lat_wgts
        1574             :     real(r8) :: to_lats(pcols), wrk(pcols)
        1575       49152 :     real(r8), allocatable, target :: wrk2d(:,:)
        1576       49152 :     real(r8), pointer :: wrk2d_in(:,:)
        1577             :     integer :: c, k, ierr, ncols
        1578             :
        1579       49152 :      nullify(wrk2d_in)
        1580      196608 :      allocate( wrk2d(cnt(1),cnt(2)), stat=ierr )
        1581       49152 :      if( ierr /= 0 ) then
        1582           0 :         write(iulog,*) 'read_2d_trc: wrk2d allocation error = ',ierr
        1583           0 :         call endrun
        1584             :      end if
        1585             :
        1586       49152 :      if(order(1)/=1 .or. order(2)/=2 .or. cnt(1)/=file%nlat .or. cnt(2)/=file%nlev) then
        1587           0 :         allocate( wrk2d_in(file%nlat, file%nlev), stat=ierr )
        1588           0 :         if( ierr /= 0 ) then
        1589           0 :            write(iulog,*) 'read_2d_trc: wrk2d_in allocation error = ',ierr
        1590           0 :            call endrun
        1591             :         end if
        1592             :      end if
        1593             :
        1594             :
        1595       49152 :     ierr = pio_get_var( fid, vid, strt, cnt, wrk2d )
        1596       49152 :     if(associated(wrk2d_in)) then
        1597           0 :        wrk2d_in = reshape( wrk2d(:,:),(/file%nlat,file%nlev/), order=order )
        1598           0 :        deallocate(wrk2d)
        1599             :     else
        1600       49152 :        wrk2d_in => wrk2d
        1601             :     end if
        1602             :
        1603      247296 :     do c=begchunk,endchunk
        1604      198144 :        ncols = get_ncols_p(c)
        1605      198144 :        call get_rlat_all_p(c, pcols, to_lats)
        1606             :
        1607      198144 :        call lininterp_init(file%lats, file%nlat, to_lats, ncols, 1, lat_wgts)
        1608    14068224 :        do k=1,file%nlev
        1609    13870080 :           call lininterp(wrk2d_in(:,k), file%nlat, wrk(1:ncols), ncols, lat_wgts)
        1610   231796224 :           loc_arr(1:ncols,k,c-begchunk+1) = wrk(1:ncols)
        1611             :        end do
        1612      247296 :        call lininterp_finish(lat_wgts)
        1613             :     end do
        1614             :
        1615       49152 :     if(allocated(wrk2d)) then
        1616       49152 :        deallocate(wrk2d)
        1617             :     else
        1618           0 :        deallocate(wrk2d_in)
        1619             :     end if
        1620       98304 :   end subroutine read_za_trc
        1621             :
        1622             : !------------------------------------------------------------------------
        1623             : ! this assumes the input data is gridded to match the physics grid
        1624           0 :   subroutine read_physgrid_2d(ncid, varname, recno, data )
        1625             :
        1626       49152 :     use ncdio_atm,        only: infld
        1627             :     use cam_grid_support, only: cam_grid_check, cam_grid_id, cam_grid_get_dim_names
        1628             :
        1629             :     type(file_desc_t) :: ncid
        1630             :     character(len=*), intent(in) :: varname
        1631             :     integer, intent(in) :: recno
        1632             :     real(r8), intent(out) :: data(1:pcols,begchunk:endchunk)
        1633             :
        1634             :     logical :: found
        1635             :     character(len=8) :: dim1name, dim2name
        1636             :     integer :: grid_id  ! grid ID for data mapping
        1637             :
        1638           0 :     grid_id = cam_grid_id('physgrid')
        1639           0 :     if (.not. cam_grid_check(grid_id)) then
        1640           0 :       call endrun('tracer_data::read_physgrid_2d: Internal error, no "physgrid" grid')
        1641             :     end if
        1642           0 :     call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
        1643             :
        1644             :     call infld( varname, ncid, dim1name, dim2name, 1, pcols, begchunk, endchunk, &
        1645           0 :                 data, found, gridname='physgrid', timelevel=recno )
        1646             :
        1647           0 :     if(.not. found) then
        1648           0 :        call endrun('tracer_data::read_physgrid_2d: Could not find '//trim(varname)//' field in input datafile')
        1649             :     end if
        1650             :
        1651           0 :   end subroutine read_physgrid_2d
        1652             :
        1653             : !------------------------------------------------------------------------
        1654             : !------------------------------------------------------------------------
        1655             : ! this assumes the input data is gridded to match the physics grid
        1656           0 :   subroutine read_physgrid_3d(ncid, varname, vrt_coord_name, nlevs, recno, data )
        1657             :
        1658           0 :     use ncdio_atm,        only: infld
        1659             :     use cam_grid_support, only: cam_grid_check, cam_grid_id, cam_grid_get_dim_names
        1660             :
        1661             :     type(file_desc_t) :: ncid
        1662             :     character(len=*), intent(in) :: varname
        1663             :     character(len=*), intent(in) :: vrt_coord_name
        1664             :     integer, intent(in) :: nlevs
        1665             :     integer, intent(in) :: recno
        1666             :     real(r8), intent(out) :: data(1:pcols,1:nlevs,begchunk:endchunk)
        1667             :
        1668             :     logical :: found
        1669             :     character(len=8) :: dim1name, dim2name
        1670             :     integer :: grid_id  ! grid ID for data mapping
        1671             :
        1672           0 :     grid_id = cam_grid_id('physgrid')
        1673           0 :     if (.not. cam_grid_check(grid_id)) then
        1674           0 :       call endrun('tracer_data::read_physgrid_3d: Internal error, no "physgrid" grid')
        1675             :     end if
        1676           0 :     call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
        1677             :
        1678             :     call infld( varname, ncid, dim1name, vrt_coord_name, dim2name, 1, pcols, 1, nlevs, begchunk, endchunk, &
        1679           0 :                 data, found, gridname='physgrid', timelevel=recno )
        1680             :
        1681           0 :     if(.not. found) then
        1682           0 :        call endrun('tracer_data::read_physgrid_3d: Could not find '//trim(varname)//' field in input datafile')
        1683             :     end if
        1684             :
        1685           0 :   end subroutine read_physgrid_3d
        1686             :
        1687             :   !------------------------------------------------------------------------
        1688             :
        1689       39936 :   subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order)
        1690           0 :     use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
        1691             :     use ppgrid,           only : pcols, begchunk, endchunk
        1692             :     use phys_grid,        only : get_ncols_p, get_rlat_all_p, get_rlon_all_p
        1693             :     use mo_constants,     only : pi
        1694             :     use dycore,           only : dycore_is
        1695             :     use polar_avg,        only : polar_average
        1696             :     use horizontal_interpolate, only : xy_interp
        1697             :
        1698             :     implicit none
        1699             :
        1700             :     type(file_desc_t), intent(in) :: fid
        1701             :     type(var_desc_t), intent(in) :: vid
        1702             :     integer, intent(in) :: strt(:), cnt(:), order(3)
        1703             :     real(r8),intent(out)  :: loc_arr(:,:,:)
        1704             :
        1705             :     type (trfile), intent(in) :: file
        1706             :
        1707             :     integer :: astat, c, ncols
        1708             :     integer :: lons(pcols), lats(pcols)
        1709             :
        1710             :     integer :: ierr
        1711             :
        1712       39936 :     real(r8), allocatable, target :: wrk3d(:,:,:)
        1713       39936 :     real(r8), pointer :: wrk3d_in(:,:,:)
        1714             :     real(r8) :: to_lons(pcols), to_lats(pcols)
        1715             :     real(r8), parameter :: zero=0_r8, twopi=2_r8*pi
        1716             :     type(interp_type) :: lon_wgts, lat_wgts
        1717             :
        1718    71359392 :     loc_arr(:,:,:) = 0._r8
        1719       39936 :     nullify(wrk3d_in)
        1720      199680 :     allocate(wrk3d(cnt(1),cnt(2),cnt(3)), stat=ierr)
        1721       39936 :     if( ierr /= 0 ) then
        1722           0 :        write(iulog,*) 'read_3d_trc: wrk3d allocation error = ',ierr
        1723           0 :        call endrun
        1724             :     end if
        1725             :
        1726       39936 :     ierr = pio_get_var( fid, vid, strt, cnt, wrk3d )
        1727             :
        1728             :     if(order(1)/=1 .or. order(2)/=2 .or. order(3)/=3 .or. &
        1729       39936 :          cnt(1)/=file%nlon.or.cnt(2)/=file%nlat.or.cnt(3)/=file%nlev) then
        1730           0 :        allocate(wrk3d_in(file%nlon,file%nlat,file%nlev),stat=ierr)
        1731           0 :        if( ierr /= 0 ) then
        1732           0 :           write(iulog,*) 'read_3d_trc: wrk3d allocation error = ',ierr
        1733           0 :           call endrun
        1734             :        end if
        1735           0 :        wrk3d_in = reshape( wrk3d(:,:,:),(/file%nlon,file%nlat,file%nlev/), order=order )
        1736           0 :        deallocate(wrk3d)
        1737             :     else
        1738       39936 :        wrk3d_in => wrk3d
        1739             :     end if
        1740             :
        1741             : ! If weighting by latitude, then perform horizontal interpolation by using weight_x, weight_y
        1742             :
        1743       39936 :    if(file%weight_by_lat) then
        1744             :
        1745           0 :      call t_startf('xy_interp')
        1746           0 :      if( file%dist ) then
        1747           0 :       do c = begchunk,endchunk
        1748           0 :         ncols = get_ncols_p(c)
        1749           0 :         lons(:ncols) = lon_global_grid_ndx(:ncols,c)
        1750           0 :         lats(:ncols) = lat_global_grid_ndx(:ncols,c)
        1751             :
        1752             :         call xy_interp(file%nlon,file%nlat,file%nlev,plon,plat,pcols,ncols, &
        1753           0 :                        file%weight0_x,file%weight0_y,wrk3d_in,loc_arr(:,:,c-begchunk+1),  &
        1754           0 :                        lons,lats,file%count0_x,file%count0_y,file%index0_x,file%index0_y)
        1755             :       enddo
        1756             :      else
        1757           0 :       do c = begchunk,endchunk
        1758           0 :         ncols = get_ncols_p(c)
        1759           0 :         lons(:ncols) = lon_global_grid_ndx(:ncols,c)
        1760           0 :         lats(:ncols) = lat_global_grid_ndx(:ncols,c)
        1761             :
        1762             :         call xy_interp(file%nlon,file%nlat,file%nlev,plon,plat,pcols,ncols,&
        1763           0 :                        file%weight_x,file%weight_y,wrk3d_in,loc_arr(:,:,c-begchunk+1), &
        1764           0 :                        lons,lats,file%count_x,file%count_y,file%index_x,file%index_y)
        1765             :       enddo
        1766             :      endif
        1767           0 :      call t_stopf('xy_interp')
        1768             :
        1769             :    else
        1770      200928 :     do c=begchunk,endchunk
        1771      160992 :        ncols = get_ncols_p(c)
        1772      160992 :        call get_rlat_all_p(c, pcols, to_lats)
        1773      160992 :        call get_rlon_all_p(c, pcols, to_lons)
        1774             :
        1775      160992 :        call lininterp_init(file%lons, file%nlon, to_lons(1:ncols), ncols, 2, lon_wgts, zero, twopi)
        1776      160992 :        call lininterp_init(file%lats, file%nlat, to_lats(1:ncols), ncols, 1, lat_wgts)
        1777             :
        1778             :
        1779      160992 :        call lininterp(wrk3d_in, file%nlon, file%nlat, file%nlev, loc_arr(:,:,c-begchunk+1), ncols, pcols, lon_wgts, lat_wgts)
        1780             :
        1781             :
        1782      160992 :        call lininterp_finish(lon_wgts)
        1783      200928 :        call lininterp_finish(lat_wgts)
        1784             :     end do
        1785             :    endif
        1786             :
        1787       39936 :     if(allocated(wrk3d)) then
        1788       39936 :        deallocate( wrk3d, stat=astat )
        1789             :     else
        1790           0 :        deallocate( wrk3d_in, stat=astat )
        1791             :     end if
        1792       39936 :     if( astat/= 0 ) then
        1793           0 :        write(iulog,*) 'read_3d_trc: failed to deallocate wrk3d array; error = ',astat
        1794           0 :        call endrun
        1795             :     endif
        1796       39936 :     if(dycore_is('LR')) call polar_average(file%nlev, loc_arr)
        1797       79872 :   end subroutine read_3d_trc
        1798             :
        1799             : !------------------------------------------------------------------------------
        1800             :
        1801     1483776 :   subroutine interpolate_trcdata( state, flds, file, pbuf2d )
        1802       39936 :     use mo_util,      only : rebin
        1803             :     use physics_types,only : physics_state
        1804             :     use physconst,    only : cday, rga
        1805             :
        1806             :     implicit none
        1807             :
        1808             :     type(physics_state), intent(in) :: state(begchunk:endchunk)
        1809             :     type (trfld),        intent(inout) :: flds(:)
        1810             :     type (trfile),       intent(inout) :: file
        1811             :
        1812             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
        1813             :
        1814             :
        1815             :     real(r8) :: fact1, fact2
        1816             :     real(r8) :: deltat
        1817             :     integer :: f,nflds,c,ncol, i,k
        1818             :     real(r8) :: ps(pcols)
        1819     2967552 :     real(r8) :: datain(pcols,file%nlev)
        1820     2967552 :     real(r8) :: pin(pcols,file%nlev)
        1821             :     real(r8)            :: model_z(pverp)
        1822             :     real(r8), parameter :: m2km  = 1.e-3_r8
        1823     1483776 :     real(r8), pointer :: data_out3d(:,:,:)
        1824     1483776 :     real(r8), pointer :: data_out(:,:)
        1825             :     integer :: chnk_offset
        1826             :     real(r8) :: data_col(pver)
        1827             :
        1828     1483776 :     nflds = size(flds)
        1829             :
        1830     1483776 :     if ( file%interp_recs == 4 ) then
        1831           0 :        deltat = file%datatimes(3) - file%datatimes(1)
        1832           0 :        fact1 = (file%datatimes(3) - file%datatimem)/deltat
        1833           0 :        fact2 = 1._r8-fact1
        1834             : !$OMP PARALLEL DO PRIVATE (C, NCOL, F)
        1835           0 :        do c = begchunk,endchunk
        1836           0 :           ncol = state(c)%ncol
        1837           0 :           if ( file%has_ps ) then
        1838           0 :              file%ps_in(1)%data(:ncol,c) = fact1*file%ps_in(1)%data(:ncol,c) + fact2*file%ps_in(3)%data(:ncol,c)
        1839             :           endif
        1840           0 :           do f = 1,nflds
        1841           0 :              flds(f)%input(1)%data(:ncol,:,c) = fact1*flds(f)%input(1)%data(:ncol,:,c) + fact2*flds(f)%input(3)%data(:ncol,:,c)
        1842             :           enddo
        1843             :        enddo
        1844             :
        1845           0 :        deltat = file%datatimes(4) - file%datatimes(2)
        1846           0 :        fact1 = (file%datatimes(4) - file%datatimep)/deltat
        1847           0 :        fact2 = 1._r8-fact1
        1848             :
        1849             : !$OMP PARALLEL DO PRIVATE (C, NCOL, F)
        1850           0 :        do c = begchunk,endchunk
        1851           0 :           ncol = state(c)%ncol
        1852           0 :           if ( file%has_ps ) then
        1853           0 :              file%ps_in(2)%data(:ncol,c) = fact1*file%ps_in(2)%data(:ncol,c) + fact2*file%ps_in(4)%data(:ncol,c)
        1854             :           endif
        1855           0 :           do f = 1,nflds
        1856           0 :              flds(f)%input(2)%data(:ncol,:,c) = fact1*flds(f)%input(2)%data(:ncol,:,c) + fact2*flds(f)%input(4)%data(:ncol,:,c)
        1857             :           enddo
        1858             :        enddo
        1859             :
        1860             :     endif
        1861             :     !-------------------------------------------------------------------------
        1862             :     !       If file%interp_recs=1 then no time interpolation -- set
        1863             :     !       fact1=1 and fact2=0 and will just use first value unmodified
        1864             :     !-------------------------------------------------------------------------
        1865             :
        1866     1483776 :     if (file%interp_recs == 1) then
        1867             :        fact1=1._r8
        1868             :        fact2=0._r8
        1869             :     else
        1870     1483776 :        file%interp_recs = 2
        1871             :
        1872     1483776 :        deltat = file%datatimep - file%datatimem
        1873             :
        1874     1483776 :        if ( file%cyclical .and. (deltat < 0._r8) ) then
        1875      741888 :           deltat = deltat+file%one_yr
        1876      741888 :           if ( file%datatimep >= file%curr_mod_time ) then
        1877      741888 :              fact1 = (file%datatimep - file%curr_mod_time)/deltat
        1878             :           else
        1879           0 :              fact1 = (file%datatimep+file%one_yr - file%curr_mod_time)/deltat
        1880             :           endif
        1881             :        else
        1882      741888 :              fact1 = (file%datatimep - file%curr_mod_time)/deltat
        1883             :        endif
        1884             :
        1885             :        ! this assures that FIXED data are b4b on restarts
        1886     1483776 :        if ( file%fixed ) then
        1887           0 :           fact1 = dble(int(fact1*cday+.5_r8))/dble(cday)
        1888             :        endif
        1889     1483776 :        fact2 = 1._r8-fact1
        1890             :     endif
        1891             :
        1892     1483776 :     chnk_offset=-begchunk+1
        1893             :
        1894    14466816 :     fld_loop: do f = 1,nflds
        1895             :
        1896    12983040 :        if (flds(f)%pbuf_ndx<=0) then
        1897     5193216 :           data_out3d => flds(f)%data(:,:,:)
        1898             :        endif
        1899             :
        1900             : !$OMP PARALLEL DO PRIVATE (C, NCOL, PS, I, K, PIN, DATAIN, MODEL_Z, DATA_OUT, DATA_COL)
        1901    66804696 :        do c = begchunk,endchunk
        1902    52337880 :           if (flds(f)%pbuf_ndx>0) then
        1903             :              call pbuf_get_field(pbuf2d, c, flds(f)%pbuf_ndx, data_out)
        1904             :           else
        1905    20935152 :              data_out => data_out3d(:,:,c+chnk_offset)
        1906             :           endif
        1907    52337880 :           ncol = state(c)%ncol
        1908    65320920 :           if (file%alt_data) then
        1909             :
        1910           0 :              if (fact2 == 0) then  ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps)
        1911           0 :                 datain(:ncol,:) = fact1*flds(f)%input(nm)%data(:ncol,:,c)
        1912             :              else
        1913           0 :                 datain(:ncol,:) = fact1*flds(f)%input(nm)%data(:ncol,:,c) + fact2*flds(f)%input(np)%data(:ncol,:,c)
        1914             :              end if
        1915           0 :              do i = 1,ncol
        1916           0 :                 model_z(1:pverp) = m2km * state(c)%zi(i,pverp:1:-1)
        1917           0 :                 if (file%geop_alt) then
        1918           0 :                    model_z(1:pverp) = model_z(1:pverp) + m2km * state(c)%phis(i)*rga
        1919             :                 endif
        1920           0 :                 if (file%conserve_column) then
        1921           0 :                    call interpz_conserve( file%nlev, pver, file%ilevs, model_z, datain(i,:), data_col(:) )
        1922             :                 else
        1923           0 :                    call rebin( file%nlev, pver, file%ilevs, model_z, datain(i,:), data_col(:) )
        1924             :                 end if
        1925           0 :                 data_out(i,:) = data_col(pver:1:-1)
        1926             :              enddo
        1927             :
        1928             :           else
        1929             :
        1930    52337880 :              if ( file%nlev>1 ) then
        1931    31402728 :                 if ( file%has_ps ) then
        1932    31402728 :                    if (fact2 == 0) then  ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps)
        1933      827136 :                       ps(:ncol) = fact1*file%ps_in(nm)%data(:ncol,c)
        1934             :                    else
        1935   523525392 :                       ps(:ncol) = fact1*file%ps_in(nm)%data(:ncol,c) + fact2*file%ps_in(np)%data(:ncol,c)
        1936             :                    end if
        1937   524352528 :                    do i = 1,ncol
        1938 21603824928 :                       do k = 1,file%nlev
        1939 21572422200 :                          pin(i,k) = file%p0*file%hyam(k) + ps(i)*file%hybm(k)
        1940             :                       enddo
        1941             :                    enddo
        1942             :                 else
        1943           0 :                    do k = 1,file%nlev
        1944           0 :                       pin(:,k) = file%levs(k)
        1945             :                    enddo
        1946             :                 endif
        1947             :              endif
        1948             :
        1949    52337880 :              if (flds(f)%srf_fld) then
        1950   349568352 :                 do i = 1,ncol
        1951   349568352 :                    if (fact2 == 0) then  ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps)
        1952           0 :                       data_out(i,1) = &
        1953           0 :                            fact1*flds(f)%input(nm)%data(i,1,c)
        1954             :                    else
        1955           0 :                       data_out(i,1) = &
        1956   328633200 :                            fact1*flds(f)%input(nm)%data(i,1,c) + fact2*flds(f)%input(np)%data(i,1,c)
        1957             :                    endif
        1958             :                 enddo
        1959             :              else
        1960    31402728 :                 if (fact2 == 0) then  ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps)
        1961    57949056 :                    datain(:ncol,:) = fact1*flds(f)%input(nm)%data(:ncol,:,c)
        1962             :                 else
        1963 22395766536 :                    datain(:ncol,:) = fact1*flds(f)%input(nm)%data(:ncol,:,c) + fact2*flds(f)%input(np)%data(:ncol,:,c)
        1964             :                 end if
        1965    31402728 :                 if ( file%top_bndry ) then
        1966           0 :                    call vert_interp_ub(ncol, file%nlev, file%levs,  datain(:ncol,:), data_out(:ncol,1) )
        1967    31402728 :                 else if ( file%top_layer ) then
        1968           0 :                    call vert_interp_ub_var(ncol, file%nlev, file%levs, state(c)%pmid(:ncol,1), datain(:ncol,:), data_out(:ncol,1) )
        1969    31402728 :                 else if(file%conserve_column) then
        1970             :                    call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, &
        1971             :                         datain, data_out(:,:), &
        1972           0 :                         file%p0,ps,file%hyai,file%hybi,file%dist)
        1973             :                 else
        1974    31402728 :                    call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, data_out(:,:) )
        1975             :                 endif
        1976             :              endif
        1977             :
        1978             :           endif
        1979             :        enddo
        1980             :
        1981             :     enddo fld_loop
        1982             :
        1983     2967552 :   end subroutine interpolate_trcdata
        1984             :
        1985             : !-----------------------------------------------------------------------
        1986             : !-----------------------------------------------------------------------
        1987       21504 :   subroutine get_dimension( fid, dname, dsize, dimid, data )
        1988             :     implicit none
        1989             :     type(file_desc_t), intent(inout) :: fid
        1990             :     character(*), intent(in) :: dname
        1991             :     integer, intent(out) :: dsize
        1992             :
        1993             :     integer, optional, intent(out) :: dimid
        1994             :     real(r8), optional, pointer, dimension(:) :: data
        1995             :
        1996             :     integer :: vid, ierr, id
        1997             :     integer :: err_handling
        1998             :
        1999       21504 :     call pio_seterrorhandling( fid, PIO_BCAST_ERROR, oldmethod=err_handling)
        2000       21504 :     ierr = pio_inq_dimid( fid, dname, id )
        2001       21504 :     call pio_seterrorhandling( fid, err_handling)
        2002             :
        2003       21504 :     if ( ierr==PIO_NOERR ) then
        2004             :
        2005       19968 :        ierr = pio_inq_dimlen( fid, id, dsize )
        2006             :
        2007       19968 :        if ( present(dimid) ) then
        2008       13824 :           dimid = id
        2009             :        endif
        2010             :
        2011       19968 :        if ( present(data) ) then
        2012       13824 :           if ( associated(data) ) then
        2013           0 :              deallocate(data, stat=ierr)
        2014             :              if( ierr /= 0 ) then
        2015             :                 write(iulog,*) 'get_dimension: data deallocation error = ',ierr
        2016             :                 call endrun('get_dimension: failed to deallocate data array')
        2017             :              end if
        2018             :           endif
        2019       41472 :           allocate( data(dsize), stat=ierr )
        2020       13824 :           if( ierr /= 0 ) then
        2021           0 :              write(iulog,*) 'get_dimension: data allocation error = ',ierr
        2022           0 :              call endrun('get_dimension: failed to allocate data array')
        2023             :           end if
        2024             :
        2025       13824 :           ierr =  pio_inq_varid( fid, dname, vid )
        2026       13824 :           ierr =  pio_get_var( fid, vid, data )
        2027             :        endif
        2028             :     else
        2029        1536 :        dsize = 1
        2030        1536 :        if ( present(dimid) ) then
        2031        1536 :           dimid = -1
        2032             :        endif
        2033             :     endif
        2034             :
        2035     1483776 :   end subroutine get_dimension
        2036             :
        2037             : !-----------------------------------------------------------------------
        2038             : !-----------------------------------------------------------------------
        2039           0 :   subroutine set_cycle_indices( fileid, cyc_ndx_beg, cyc_ndx_end, cyc_yr )
        2040             :
        2041             :     implicit none
        2042             :
        2043             :     type(file_desc_t), intent(inout)  :: fileid
        2044             :     integer, intent(out) :: cyc_ndx_beg
        2045             :     integer, intent(out) :: cyc_ndx_end
        2046             :     integer, intent(in)  :: cyc_yr
        2047             :
        2048           0 :     integer, allocatable , dimension(:) :: dates
        2049             :     integer :: timesize, i, astat, year, ierr
        2050             :     type(var_desc_T) :: dateid
        2051           0 :     call get_dimension( fileid, 'time', timesize )
        2052           0 :     cyc_ndx_beg=-1
        2053             :
        2054           0 :     allocate( dates(timesize), stat=astat  )
        2055           0 :     if( astat/= 0 ) then
        2056           0 :        write(*,*) 'set_cycle_indices: failed to allocate dates array; error = ',astat
        2057           0 :        call endrun
        2058             :     end if
        2059             :
        2060           0 :     ierr = pio_inq_varid(   fileid, 'date',  dateid  )
        2061           0 :     ierr = pio_get_var( fileid, dateid, dates )
        2062             :
        2063           0 :     do i=1,timesize
        2064           0 :        year = dates(i) / 10000
        2065           0 :        if ( year == cyc_yr ) then
        2066           0 :           if  (cyc_ndx_beg < 0)  then
        2067           0 :              cyc_ndx_beg = i
        2068             :           endif
        2069           0 :           cyc_ndx_end = i
        2070             :        endif
        2071             :     enddo
        2072           0 :     deallocate( dates, stat=astat  )
        2073           0 :     if( astat/= 0 ) then
        2074           0 :        write(*,*) 'set_cycle_indices: failed to deallocate dates array; error = ',astat
        2075           0 :        call endrun
        2076             :     end if
        2077           0 :     if (cyc_ndx_beg < 0) then
        2078           0 :        write(*,*) 'set_cycle_indices: cycle year not found : ' , cyc_yr
        2079           0 :        call endrun('set_cycle_indices: cycle year not found')
        2080             :     endif
        2081             :
        2082           0 :   end subroutine set_cycle_indices
        2083             : !-----------------------------------------------------------------------
        2084             :
        2085             : !-----------------------------------------------------------------------
        2086        6144 :   subroutine open_trc_datafile( fname, path, piofile, times, cyc_ndx_beg, cyc_ndx_end, cyc_yr )
        2087             :
        2088             :     use ioFileMod,     only: getfil
        2089             :     use cam_pio_utils, only: cam_pio_openfile
        2090             :
        2091             :     implicit none
        2092             :
        2093             :     character(*), intent(in) :: fname
        2094             :     character(*), intent(in) :: path
        2095             :     type(file_desc_t), intent(inout) :: piofile
        2096             :     real(r8), pointer :: times(:)
        2097             :
        2098             :     integer, optional, intent(out) :: cyc_ndx_beg
        2099             :     integer, optional, intent(out) :: cyc_ndx_end
        2100             :     integer, optional, intent(in) :: cyc_yr
        2101             :
        2102             :     character(len=shr_kind_cl) :: filen, filepath
        2103             :     integer :: year, month, day, i, timesize
        2104             :     integer :: dateid,secid
        2105        6144 :     integer, allocatable , dimension(:) :: dates, datesecs
        2106             :     integer :: astat, ierr
        2107             :     logical :: need_first_ndx
        2108             :     integer :: err_handling
        2109             :
        2110        6144 :     if (len_trim(path) == 0) then
        2111           0 :        filepath = trim(fname)
        2112             :     else
        2113        6144 :        filepath = trim(path) // '/' // trim(fname)
        2114             :     end if
        2115             :     !
        2116             :     ! open file and get fileid
        2117             :     !
        2118        6144 :     call getfil( filepath, filen, 0 )
        2119        6144 :     call cam_pio_openfile( piofile, filen, PIO_NOWRITE)
        2120        6144 :     if(masterproc) write(iulog,*)'open_trc_datafile: ',trim(filen)
        2121             :
        2122        6144 :     call get_dimension(piofile, 'time', timesize)
        2123             :
        2124        6144 :     if ( associated(times) ) then
        2125           0 :        deallocate(times, stat=ierr)
        2126             :        if( ierr /= 0 ) then
        2127             :           write(iulog,*) 'open_trc_datafile: data deallocation error = ',ierr
        2128             :           call endrun('open_trc_datafile: failed to deallocate data array')
        2129             :        end if
        2130             :     endif
        2131       18432 :     allocate( times(timesize), stat=ierr )
        2132        6144 :     if( ierr /= 0 ) then
        2133           0 :        write(iulog,*) 'open_trc_datafile: data allocation error = ',ierr
        2134           0 :        call endrun('open_trc_datafile: failed to allocate data array')
        2135             :     end if
        2136             :
        2137       18432 :     allocate( dates(timesize), stat=astat  )
        2138        6144 :     if( astat/= 0 ) then
        2139           0 :        if(masterproc) write(iulog,*) 'open_trc_datafile: failed to allocate dates array; error = ',astat
        2140           0 :        call endrun
        2141             :     end if
        2142       18432 :     allocate( datesecs(timesize), stat=astat  )
        2143        6144 :     if( astat/= 0 ) then
        2144           0 :        if(masterproc) write(iulog,*) 'open_trc_datafile: failed to allocate datesec array; error = ',astat
        2145           0 :        call endrun
        2146             :     end if
        2147             :
        2148        6144 :     ierr =  pio_inq_varid( piofile, 'date',    dateid  )
        2149        6144 :     call pio_seterrorhandling( piofile, PIO_BCAST_ERROR, oldmethod=err_handling)
        2150        6144 :     ierr = pio_inq_varid( piofile, 'datesec', secid  )
        2151        6144 :     call pio_seterrorhandling( piofile, err_handling)
        2152             :
        2153        6144 :     if(ierr==PIO_NOERR) then
        2154        4608 :        ierr = pio_get_var( piofile, secid,  datesecs  )
        2155             :     else
        2156     2913792 :        datesecs=0
        2157             :     end if
        2158             :
        2159        6144 :     ierr =  pio_get_var( piofile, dateid, dates )
        2160        6144 :     need_first_ndx=.true.
        2161             :
        2162    39945216 :     do i=1,timesize
        2163    39939072 :        year = dates(i) / 10000
        2164    39939072 :        month = mod(dates(i),10000)/100
        2165    39939072 :        day = mod(dates(i),100)
        2166    39939072 :        call set_time_float_from_date( times(i), year, month, day, datesecs(i) )
        2167    39945216 :        if ( present(cyc_yr) ) then
        2168     2930688 :           if ( year == cyc_yr ) then
        2169       36864 :              if ( present(cyc_ndx_beg) .and. need_first_ndx ) then
        2170        3072 :                 cyc_ndx_beg = i
        2171        3072 :                 need_first_ndx = .false.
        2172             :              endif
        2173       36864 :              if ( present(cyc_ndx_end) ) then
        2174       36864 :                 cyc_ndx_end = i
        2175             :              endif
        2176             :           endif
        2177             :        endif
        2178             :     enddo
        2179             :
        2180        6144 :     deallocate( dates, stat=astat  )
        2181        6144 :     if( astat/= 0 ) then
        2182           0 :        if(masterproc) write(iulog,*) 'open_trc_datafile: failed to deallocate dates array; error = ',astat
        2183           0 :        call endrun
        2184             :     end if
        2185        6144 :     deallocate( datesecs, stat=astat  )
        2186        6144 :     if( astat/= 0 ) then
        2187           0 :        if(masterproc) write(iulog,*) 'open_trc_datafile: failed to deallocate datesec array; error = ',astat
        2188           0 :        call endrun
        2189             :     end if
        2190             :
        2191        6144 :     if ( present(cyc_yr) .and. present(cyc_ndx_beg) ) then
        2192        3072 :        if (cyc_ndx_beg < 0) then
        2193           0 :           write(iulog,*) 'open_trc_datafile: cycle year not found : ' , cyc_yr
        2194           0 :           call endrun('open_trc_datafile: cycle year not found '//trim(filepath))
        2195             :        endif
        2196             :     endif
        2197             :
        2198       12288 :   end subroutine open_trc_datafile
        2199             :
        2200             : !--------------------------------------------------------------------------
        2201             : !--------------------------------------------------------------------------
        2202        6144 :   subroutine specify_fields( specifier, fields )
        2203             :
        2204             :     implicit none
        2205             :
        2206             :     character(len=*), intent(in) :: specifier(:)
        2207             :     type(trfld), pointer, dimension(:) :: fields
        2208             :
        2209             :     integer :: fld_cnt, astat
        2210             :     integer :: i,j
        2211             :     character(len=shr_kind_cl) :: str1, str2
        2212        6144 :     character(len=32), allocatable, dimension(:) :: fld_name,  src_name
        2213             :     integer :: nflds
        2214             :
        2215        6144 :     nflds = size(specifier)
        2216             :
        2217       24576 :     allocate(fld_name(nflds),  src_name(nflds), stat=astat )
        2218        6144 :     if( astat/= 0 ) then
        2219           0 :        write(iulog,*) 'specify_fields: failed to allocate fld_name, src_name arrays; error = ',astat
        2220           0 :        call endrun
        2221             :     end if
        2222             :
        2223             :     fld_cnt = 0
        2224             :
        2225       59904 :     count_cnst: do i = 1, nflds
        2226             :
        2227       56832 :        if ( len_trim( specifier(i) ) == 0 ) then
        2228             :           exit count_cnst
        2229             :        endif
        2230             :
        2231       53760 :        j = scan( specifier(i),':')
        2232             :
        2233       53760 :        if (j > 0) then
        2234       32256 :           str1 = trim(adjustl( specifier(i)(:j-1) ))
        2235       32256 :           str2 = trim(adjustl( specifier(i)(j+1:) ))
        2236       32256 :           fld_name(i) = trim(adjustl( str1 ))
        2237       32256 :           src_name(i) = trim(adjustl( str2 ))
        2238             :        else
        2239       21504 :           fld_name(i) = trim(adjustl( specifier(i) ))
        2240       21504 :           src_name(i) = trim(adjustl( specifier(i) ))
        2241             :        endif
        2242             :
        2243       59904 :        fld_cnt = fld_cnt + 1
        2244             :
        2245             :     enddo count_cnst
        2246             :
        2247        6144 :     if( fld_cnt < 1 ) then
        2248           0 :        nullify(fields)
        2249           0 :        return
        2250             :     end if
        2251             :
        2252             :     !-----------------------------------------------------------------------
        2253             :     !   ... allocate field type array
        2254             :     !-----------------------------------------------------------------------
        2255       96768 :     allocate( fields(fld_cnt), stat=astat )
        2256        6144 :     if( astat/= 0 ) then
        2257           0 :        write(iulog,*) 'specify_fields: failed to allocate fields array; error = ',astat
        2258           0 :        call endrun
        2259             :     end if
        2260             :
        2261       59904 :     do i = 1,fld_cnt
        2262       53760 :        fields(i)%fldnam = fld_name(i)
        2263       59904 :        fields(i)%srcnam = src_name(i)
        2264             :     enddo
        2265             :
        2266        6144 :     deallocate(fld_name, src_name)
        2267             :
        2268        6144 :   end subroutine specify_fields
        2269             :
        2270             : !------------------------------------------------------------------------------
        2271             :
        2272        6144 :   subroutine init_trc_restart( whence, piofile, tr_file )
        2273             :
        2274             :     implicit none
        2275             :     character(len=*), intent(in) :: whence
        2276             :     type(file_desc_t), intent(inout) :: piofile
        2277             :     type(trfile), intent(inout) :: tr_file
        2278             :
        2279             :     character(len=32) :: name
        2280             :     integer :: ioerr, mcdimid, maxlen
        2281             :     integer :: err_handling
        2282             :
        2283             :     ! Dimension should already be defined in restart file
        2284        6144 :     call pio_seterrorhandling(pioFile, PIO_BCAST_ERROR, oldmethod=err_handling)
        2285        6144 :     ioerr = pio_inq_dimid(pioFile,'max_chars', mcdimid)
        2286        6144 :     call pio_seterrorhandling(pioFile, err_handling)
        2287             :     ! but define it if nessasary
        2288        6144 :     if(ioerr/= PIO_NOERR) then
        2289        1536 :        ioerr = pio_def_dim(pioFile, 'max_chars', SHR_KIND_CL, mcdimid)
        2290             :     end if
        2291             :
        2292        6144 :     if(len_trim(tr_file%curr_filename)>1) then
        2293        3072 :        allocate(tr_file%currfnameid)
        2294        3072 :        name = trim(whence)//'_curr_fname'
        2295        6144 :        ioerr = pio_def_var(pioFile, name,pio_char,  (/mcdimid/), tr_file%currfnameid)
        2296        3072 :        ioerr = pio_put_att(pioFile, tr_file%currfnameid, 'offset_time', tr_file%offset_time)
        2297        3072 :        maxlen = len_trim(tr_file%curr_filename)
        2298        3072 :        ioerr = pio_put_att(pioFile, tr_file%currfnameid, 'actual_len', maxlen)
        2299             :     else
        2300        3072 :        nullify(tr_file%currfnameid)
        2301             :     end if
        2302             :
        2303        6144 :     if(len_trim(tr_file%next_filename)>1) then
        2304           0 :        allocate(tr_file%nextfnameid)
        2305           0 :        name = trim(whence)//'_next_fname'
        2306           0 :        ioerr = pio_def_var(pioFile, name,pio_char,  (/mcdimid/), tr_file%nextfnameid)
        2307           0 :        maxlen = len_trim(tr_file%next_filename)
        2308           0 :        ioerr = pio_put_att(pioFile, tr_file%nextfnameid, 'actual_len', maxlen)
        2309             :     else
        2310        6144 :        nullify(tr_file%nextfnameid)
        2311             :     end if
        2312        6144 :   end subroutine init_trc_restart
        2313             : !-------------------------------------------------------------------------
        2314             : ! writes file names to restart file
        2315             : !-------------------------------------------------------------------------
        2316        6144 :   subroutine write_trc_restart( piofile, tr_file )
        2317             :
        2318             :     implicit none
        2319             :
        2320             :     type(file_desc_t), intent(inout) :: piofile
        2321             :     type(trfile), intent(inout) :: tr_file
        2322             :
        2323             :     integer :: ioerr   ! error status
        2324        6144 :     if(associated(tr_file%currfnameid)) then
        2325        3072 :        ioerr = pio_put_var(pioFile, tr_file%currfnameid, tr_file%curr_filename)
        2326        3072 :        deallocate(tr_file%currfnameid)
        2327             :        nullify(tr_file%currfnameid)
        2328             :     end if
        2329        6144 :     if(associated(tr_file%nextfnameid)) then
        2330           0 :        ioerr = pio_put_var(pioFile, tr_file%nextfnameid, tr_file%next_filename)
        2331           0 :        deallocate(tr_file%nextfnameid)
        2332             :        nullify(tr_file%nextfnameid)
        2333             :     end if
        2334        6144 :   end subroutine write_trc_restart
        2335             :
        2336             : !-------------------------------------------------------------------------
        2337             : ! reads file names from restart file
        2338             : !-------------------------------------------------------------------------
        2339        3072 :   subroutine read_trc_restart( whence, piofile, tr_file )
        2340             :
        2341             :     implicit none
        2342             :
        2343             :     character(len=*), intent(in) :: whence
        2344             :     type(file_desc_t), intent(inout) :: piofile
        2345             :     type(trfile), intent(inout) :: tr_file
        2346             :     type(var_desc_t) :: vdesc
        2347             :     character(len=64) :: name
        2348             :     integer :: ioerr   ! error status
        2349             :     integer :: slen
        2350             :     integer :: err_handling
        2351             :
        2352        3072 :     call PIO_SetErrorHandling(piofile, PIO_BCAST_ERROR, oldmethod=err_handling)
        2353        3072 :     name = trim(whence)//'_curr_fname'
        2354        3072 :     ioerr = pio_inq_varid(piofile, name, vdesc)
        2355        3072 :     if(ioerr==PIO_NOERR) then
        2356        1536 :        tr_file%curr_filename=' '
        2357        1536 :        ioerr = pio_get_att(piofile, vdesc, 'offset_time', tr_file%offset_time)
        2358        1536 :        ioerr = pio_get_att(piofile, vdesc, 'actual_len', slen)
        2359        1536 :        ioerr = pio_get_var(piofile, vdesc, tr_file%curr_filename)
        2360        1536 :        if(slen<SHR_KIND_CL) tr_file%curr_filename(slen+1:)=' '
        2361             :     end if
        2362             :
        2363        3072 :     name = trim(whence)//'_next_fname'
        2364        3072 :     ioerr = pio_inq_varid(piofile, name, vdesc)
        2365        3072 :     if(ioerr==PIO_NOERR) then
        2366           0 :        tr_file%next_filename=' '
        2367           0 :        ioerr = pio_get_att(piofile, vdesc, 'actual_len', slen)
        2368           0 :        ioerr = pio_get_var(piofile, vdesc, tr_file%next_filename)
        2369           0 :        if(slen<SHR_KIND_CL) tr_file%next_filename(slen+1:)=' '
        2370             :     end if
        2371        3072 :     call PIO_SetErrorHandling(piofile, err_handling)
        2372             :
        2373             :
        2374             :
        2375        3072 :   end subroutine read_trc_restart
        2376             : !------------------------------------------------------------------------------
        2377           0 :   subroutine interpz_conserve( nsrc, ntrg, src_x, trg_x, src, trg)
        2378             :
        2379             :     implicit none
        2380             :
        2381             :     integer, intent(in)   :: nsrc                  ! dimension source array
        2382             :     integer, intent(in)   :: ntrg                  ! dimension target array
        2383             :     real(r8), intent(in)  :: src_x(nsrc+1)         ! source coordinates
        2384             :     real(r8), intent(in)  :: trg_x(ntrg+1)         ! target coordinates
        2385             :     real(r8), intent(in)  :: src(nsrc)             ! source array
        2386             :     real(r8), intent(out) :: trg(ntrg)             ! target array
        2387             :
        2388             :     !---------------------------------------------------------------
        2389             :     !   ... local variables
        2390             :     !---------------------------------------------------------------
        2391             :     integer  :: i, j
        2392             :     integer  :: sil
        2393             :     real(r8) :: tl, y
        2394             :     real(r8) :: bot, top
        2395             :
        2396             :
        2397           0 :     do i = 1, ntrg
        2398           0 :        tl = trg_x(i)
        2399           0 :        if ( (tl<src_x(nsrc+1)).and.(trg_x(i+1)>src_x(1)) ) then
        2400           0 :           do sil = 1,nsrc
        2401           0 :              if ( (tl-src_x(sil))*(tl-src_x(sil+1))<=0.0_r8 ) then
        2402             :                 exit
        2403             :              end if
        2404             :           end do
        2405             :
        2406           0 :           if ( tl<src_x(1) ) sil = 1
        2407             :
        2408           0 :           y = 0.0_r8
        2409           0 :           bot = max(tl,src_x(1))
        2410             :           top = trg_x(i+1)
        2411           0 :           do j = sil, nsrc
        2412           0 :              if ( top>src_x(j+1) ) then
        2413           0 :                 y = y+(src_x(j+1)-bot)*src(j)/(src_x(j+1)-src_x(j))
        2414             :                 bot = src_x(j+1)
        2415             :              else
        2416           0 :                 y = y+(top-bot)*src(j)/(src_x(j+1)-src_x(j))
        2417           0 :                 exit
        2418             :              endif
        2419             :           enddo
        2420           0 :           trg(i) = y
        2421             :        else
        2422           0 :           trg(i) = 0.0_r8
        2423             :        end if
        2424             :     end do
        2425             :
        2426           0 :     if ( trg_x(1)>src_x(1) ) then
        2427             :        top = trg_x(1)
        2428             :        bot = src_x(1)
        2429             :        y = 0.0_r8
        2430           0 :        do j = 1, nsrc
        2431           0 :           if ( top>src_x(j+1) ) then
        2432           0 :              y = y+(src_x(j+1)-bot)*src(j)/(src_x(j+1)-src_x(j))
        2433             :              bot = src_x(j+1)
        2434             :           else
        2435           0 :              y = y+(top-bot)*src(j)/(src_x(j+1)-src_x(j))
        2436           0 :              exit
        2437             :           endif
        2438             :        enddo
        2439           0 :        trg(1) = trg(1)+y
        2440             :     endif
        2441             :
        2442             :
        2443           0 :   end subroutine interpz_conserve
        2444             :
        2445             : !------------------------------------------------------------------------------
        2446           0 :    subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi, use_flight_distance)
        2447             :
        2448             :     implicit none
        2449             :
        2450             :     integer, intent(in)   :: ncol
        2451             :     integer, intent(in)   :: nsrc                  ! dimension source array
        2452             :     integer, intent(in)   :: ntrg                  ! dimension target array
        2453           0 :     real(r8)              :: src_x(nsrc+1)         ! source coordinates
        2454             :     real(r8), intent(in)      :: trg_x(pcols,ntrg+1)         ! target coordinates
        2455             :     real(r8), intent(in)      :: src(pcols,nsrc)             ! source array
        2456             :     logical, intent(in)   :: use_flight_distance                    ! .true. = flight distance, .false. = mixing ratio
        2457             :     real(r8), intent(out)     :: trg(pcols,ntrg)             ! target array
        2458             :
        2459             :     real(r8) :: ps(pcols), p0, hyai(nsrc+1), hybi(nsrc+1)
        2460             :     !---------------------------------------------------------------
        2461             :     !   ... local variables
        2462             :     !---------------------------------------------------------------
        2463             :     integer  :: i, j, n
        2464             :     integer  :: sil
        2465             :     real(r8)     :: tl, y
        2466             :     real(r8)     :: bot, top
        2467             :
        2468             :
        2469             :
        2470           0 :     do n = 1,ncol
        2471             :
        2472           0 :        do i=1,nsrc+1
        2473           0 :           src_x(i) = p0*hyai(i)+ps(n)*hybi(i)
        2474             :        enddo
        2475             :
        2476           0 :        do i = 1, ntrg
        2477           0 :           tl = trg_x(n,i+1)
        2478           0 :           if( (tl>src_x(1)).and.(trg_x(n,i)<src_x(nsrc+1)) ) then
        2479           0 :              do sil = 1,nsrc
        2480           0 :                 if( (tl-src_x(sil))*(tl-src_x(sil+1))<=0.0_r8 ) then
        2481             :                    exit
        2482             :                 end if
        2483             :              end do
        2484             :
        2485           0 :           if( tl>src_x(nsrc+1)) sil = nsrc
        2486             :
        2487           0 :              y = 0.0_r8
        2488           0 :              bot = min(tl,src_x(nsrc+1))
        2489             :              top = trg_x(n,i)
        2490           0 :              do j = sil,1,-1
        2491           0 :                 if( top<src_x(j) ) then
        2492           0 :                     if(use_flight_distance) then
        2493           0 :                         y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j))
        2494             :                     else
        2495           0 :                         y = y+(bot-src_x(j))*src(n,j)
        2496             :                     endif
        2497           0 :                     bot = src_x(j)
        2498             :                 else
        2499           0 :                    if(use_flight_distance) then
        2500           0 :                       y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j))
        2501             :                    else
        2502           0 :                       y = y+(bot-top)*src(n,j)
        2503             :                    endif
        2504             :                    exit
        2505             :                 endif
        2506             :              enddo
        2507           0 :              trg(n,i) = y
        2508             :           else
        2509           0 :              trg(n,i) = 0.0_r8
        2510             :           end if
        2511             :        end do
        2512             :
        2513           0 :        if( trg_x(n,ntrg+1)<src_x(nsrc+1) ) then
        2514             :           top = trg_x(n,ntrg+1)
        2515             :           bot = src_x(nsrc+1)
        2516             :           y = 0.0_r8
        2517           0 :           do j=nsrc,1,-1
        2518           0 :              if( top<src_x(j) ) then
        2519           0 :                 if(use_flight_distance) then
        2520           0 :                    y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j))
        2521             :                 else
        2522           0 :                    y = y+(bot-src_x(j))*src(n,j)
        2523             :                 endif
        2524           0 :                 bot = src_x(j)
        2525             :              else
        2526           0 :                 if(use_flight_distance) then
        2527           0 :                    y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j))
        2528             :                 else
        2529           0 :                    y = y+(bot-top)*src(n,j)
        2530             :                 endif
        2531             :                 exit
        2532             :              endif
        2533             :           enddo
        2534           0 :           trg(n,ntrg) = trg(n,ntrg)+y
        2535             :        endif
        2536             :
        2537             :        ! turn mass into mixing ratio
        2538           0 :        if(.not. use_flight_distance) then
        2539           0 :           do i=1,ntrg
        2540           0 :              trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i))
        2541             :           enddo
        2542             :        endif
        2543             :
        2544             :     enddo
        2545             :
        2546           0 :    end subroutine vert_interp_mixrat
        2547             : !------------------------------------------------------------------------------
        2548    31402728 :   subroutine vert_interp( ncol, levsiz, pin, pmid, datain, dataout )
        2549             :     !--------------------------------------------------------------------------
        2550             :     !
        2551             :     ! Interpolate data from current time-interpolated values to model levels
        2552             :     !--------------------------------------------------------------------------
        2553             :     implicit none
        2554             :     ! Arguments
        2555             :     !
        2556             :     integer,  intent(in)  :: ncol                ! number of atmospheric columns
        2557             :     integer,  intent(in)  :: levsiz
        2558             :     real(r8), intent(in)  :: pin(pcols,levsiz)
        2559             :     real(r8), intent(in)  :: pmid(pcols,pver)          ! level pressures
        2560             :     real(r8), intent(in)  :: datain(pcols,levsiz)
        2561             :     real(r8), intent(out) :: dataout(pcols,pver)
        2562             :
        2563             :     !
        2564             :     ! local storage
        2565             :     !
        2566             :
        2567             :     integer ::  i                   ! longitude index
        2568             :     integer ::  k, kk, kkstart      ! level indices
        2569             :     integer ::  kupper(pcols)       ! Level indices for interpolation
        2570             :     real(r8) :: dpu                ! upper level pressure difference
        2571             :     real(r8) :: dpl                ! lower level pressure difference
        2572             :
        2573             :
        2574             :
        2575             :     !--------------------------------------------------------------------------
        2576             :     !
        2577             :     ! Initialize index array
        2578             :     !
        2579   524352528 :     do i=1,ncol
        2580   524352528 :        kupper(i) = 1
        2581             :     end do
        2582             :
        2583  2951856432 :     do k=1,pver
        2584             :        !
        2585             :        ! Top level we need to start looking is the top level for the previous k
        2586             :        ! for all column points
        2587             :        !
        2588             :        kkstart = levsiz
        2589 48764785104 :        do i=1,ncol
        2590 48764785104 :           kkstart = min0(kkstart,kupper(i))
        2591             :        end do
        2592             :        !
        2593             :        ! Store level indices for interpolation
        2594             :        !
        2595 56832625862 :        do kk=kkstart,levsiz-1
        2596 >90314*10^7 :           do i=1,ncol
        2597 >90022*10^7 :              if (pin(i,kk)<pmid(i,k) .and. pmid(i,k)<=pin(i,kk+1)) then
        2598 41220335157 :                 kupper(i) = kk
        2599             :              end if
        2600             :           end do
        2601             :        end do
        2602             :        ! interpolate or extrapolate...
        2603 48796187832 :        do i=1,ncol
        2604 48764785104 :           if (pmid(i,k) < pin(i,1)) then
        2605  3661912800 :              dataout(i,k) = datain(i,1)*pmid(i,k)/pin(i,1)
        2606 42182418600 :           else if (pmid(i,k) > pin(i,levsiz)) then
        2607   962083443 :              dataout(i,k) = datain(i,levsiz)
        2608             :           else
        2609 41220335157 :              dpu = pmid(i,k) - pin(i,kupper(i))
        2610 41220335157 :              dpl = pin(i,kupper(i)+1) - pmid(i,k)
        2611             :              dataout(i,k) = (datain(i,kupper(i) )*dpl + &
        2612 41220335157 :                   datain(i,kupper(i)+1)*dpu)/(dpl + dpu)
        2613             :           end if
        2614             :        end do
        2615             :     end do
        2616             :
        2617             :
        2618    31402728 :   end subroutine vert_interp
        2619             :
        2620             : !------------------------------------------------------------------------------
        2621           0 :   subroutine vert_interp_ub( ncol, nlevs, plevs,  datain, dataout )
        2622             :     use ref_pres, only : ptop_ref
        2623             :
        2624             :
        2625             :     !-----------------------------------------------------------------------
        2626             :     !
        2627             :     ! Interpolate data from current time-interpolated values to top interface pressure
        2628             :     !  -- from mo_tgcm_ubc.F90
        2629             :     !--------------------------------------------------------------------------
        2630             :     implicit none
        2631             :     ! Arguments
        2632             :     !
        2633             :     integer,  intent(in)  :: ncol
        2634             :     integer,  intent(in)  :: nlevs
        2635             :     real(r8), intent(in)  :: plevs(nlevs)
        2636             :     real(r8), intent(in)  :: datain(ncol,nlevs)
        2637             :     real(r8), intent(out) :: dataout(ncol)
        2638             :
        2639             :     !
        2640             :     ! local variables
        2641             :     !
        2642             :     integer  :: i,ku,kl,kk
        2643             :     real(r8) :: pinterp, delp
        2644             :
        2645           0 :     pinterp = ptop_ref
        2646             :
        2647           0 :     if( pinterp <= plevs(1) ) then
        2648             :        kl = 1
        2649             :        ku = 1
        2650             :        delp = 0._r8
        2651           0 :     else if( pinterp >= plevs(nlevs) ) then
        2652             :        kl = nlevs
        2653             :        ku = nlevs
        2654             :        delp = 0._r8
        2655             :     else
        2656             :
        2657           0 :        do kk = 2,nlevs
        2658           0 :           if( pinterp <= plevs(kk) ) then
        2659           0 :              ku = kk
        2660           0 :              kl = kk - 1
        2661           0 :              delp = log( pinterp/plevs(kk) ) / log( plevs(kk-1)/plevs(kk) )
        2662           0 :              exit
        2663             :           end if
        2664             :        end do
        2665             :
        2666             :     end if
        2667             :
        2668           0 :     do i = 1,ncol
        2669           0 :        dataout(i) = datain(i,kl) + delp * (datain(i,ku) - datain(i,kl))
        2670             :     end do
        2671             :
        2672           0 :   end subroutine vert_interp_ub
        2673             : !------------------------------------------------------------------------------
        2674             : !------------------------------------------------------------------------------
        2675           0 :   subroutine vert_interp_ub_var( ncol, nlevs, plevs, press, datain, dataout )
        2676             :
        2677             :     !-----------------------------------------------------------------------
        2678             :     !
        2679             :     ! Interpolate data from current time-interpolated values to press
        2680             :     !
        2681             :     !--------------------------------------------------------------------------
        2682             :     ! Arguments
        2683             :     !
        2684             :     integer,  intent(in)  :: ncol
        2685             :     integer,  intent(in)  :: nlevs
        2686             :     real(r8), intent(in)  :: plevs(nlevs)
        2687             :     real(r8), intent(in)  :: press(ncol)
        2688             :     real(r8), intent(in)  :: datain(ncol,nlevs)
        2689             :     real(r8), intent(out) :: dataout(ncol)
        2690             :
        2691             :     !
        2692             :     ! local variables
        2693             :     !
        2694             :     integer  :: i,k
        2695             :     integer  :: ku,kl
        2696             :     real(r8) :: delp
        2697             :
        2698             :
        2699           0 :     do i = 1,ncol
        2700             :
        2701           0 :        if( press(i) <= plevs(1) ) then
        2702             :           kl = 1
        2703             :           ku = 1
        2704             :           delp = 0._r8
        2705           0 :        else if( press(i) >= plevs(nlevs) ) then
        2706             :           kl = nlevs
        2707             :           ku = nlevs
        2708             :           delp = 0._r8
        2709             :        else
        2710             :
        2711           0 :           do k = 2,nlevs
        2712           0 :              if( press(i) <= plevs(k) ) then
        2713           0 :                 ku = k
        2714           0 :                 kl = k - 1
        2715           0 :                 delp = log( press(i)/plevs(k) ) / log( plevs(k-1)/plevs(k) )
        2716           0 :                 exit
        2717             :              end if
        2718             :           end do
        2719             :
        2720             :        end if
        2721             :
        2722           0 :        dataout(i) = datain(i,kl) + delp * (datain(i,ku) - datain(i,kl))
        2723             :     end do
        2724             :
        2725           0 :   end subroutine vert_interp_ub_var
        2726             : !------------------------------------------------------------------------------
        2727             :
        2728             : !------------------------------------------------------------------------------
        2729             : !------------------------------------------------------------------------------
        2730           0 :   subroutine advance_file(file)
        2731             :
        2732             :     !------------------------------------------------------------------------------
        2733             :     !   This routine advances to the next file
        2734             :     !------------------------------------------------------------------------------
        2735             :
        2736             :     use shr_sys_mod, only: shr_sys_system
        2737             :     use ioFileMod, only: getfil
        2738             :
        2739             :     implicit none
        2740             :
        2741             :     type(trfile), intent(inout) :: file
        2742             :
        2743             :     !-----------------------------------------------------------------------
        2744             :     !   local variables
        2745             :     !-----------------------------------------------------------------------
        2746             :     character(len=shr_kind_cl) :: ctmp
        2747             :     character(len=shr_kind_cl) :: loc_fname
        2748             :     integer            :: istat, astat
        2749             :
        2750             :     !-----------------------------------------------------------------------
        2751             :     !   close current file ...
        2752             :     !-----------------------------------------------------------------------
        2753           0 :     call pio_closefile( file%curr_fileid )
        2754             :
        2755             :     !-----------------------------------------------------------------------
        2756             :     !   remove if requested
        2757             :     !-----------------------------------------------------------------------
        2758           0 :     if( file%remove_trc_file ) then
        2759           0 :        call getfil( file%curr_filename, loc_fname, 0 )
        2760           0 :        write(iulog,*) 'advance_file: removing file = ',trim(loc_fname)
        2761           0 :        ctmp = 'rm -f ' // trim(loc_fname)
        2762           0 :        write(iulog,*) 'advance_file: fsystem issuing command - '
        2763           0 :        write(iulog,*) trim(ctmp)
        2764           0 :        call shr_sys_system( ctmp, istat )
        2765             :     end if
        2766             :
        2767             :     !-----------------------------------------------------------------------
        2768             :     !   Advance the filename and file id
        2769             :     !-----------------------------------------------------------------------
        2770           0 :     file%curr_filename = file%next_filename
        2771           0 :     file%curr_fileid = file%next_fileid
        2772             :
        2773             :     !-----------------------------------------------------------------------
        2774             :     !   Advance the curr_data_times
        2775             :     !-----------------------------------------------------------------------
        2776           0 :     deallocate( file%curr_data_times, stat=astat )
        2777           0 :     if( astat/= 0 ) then
        2778           0 :        write(iulog,*) 'advance_file: failed to deallocate file%curr_data_times array; error = ',astat
        2779           0 :        call endrun
        2780             :     end if
        2781           0 :     allocate( file%curr_data_times( size( file%next_data_times ) ), stat=astat )
        2782           0 :     if( astat/= 0 ) then
        2783           0 :        write(iulog,*) 'advance_file: failed to allocate file%curr_data_times array; error = ',astat
        2784           0 :        call endrun
        2785             :     end if
        2786           0 :     file%curr_data_times(:) = file%next_data_times(:)
        2787             :
        2788             :     !-----------------------------------------------------------------------
        2789             :     !   delete information about next file (as was just assigned to current)
        2790             :     !-----------------------------------------------------------------------
        2791           0 :     file%next_filename = ''
        2792             :
        2793           0 :     deallocate( file%next_data_times, stat=astat )
        2794           0 :     if( astat/= 0 ) then
        2795           0 :        write(iulog,*) 'advance_file: failed to deallocate file%next_data_times array; error = ',astat
        2796           0 :        call endrun
        2797             :     end if
        2798           0 :     nullify( file%next_data_times )
        2799             :
        2800           0 :   end subroutine advance_file
        2801             :
        2802           0 : end module tracer_data

Generated by: LCOV version 1.14