LCOV - code coverage report
Current view: | top level - control - sat_hist.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 31 | 496 | 6.2 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 17 | 17.6 %

              Line data    Source code
           1             : !-------------------------------------------------------------------------------
           2             : ! Outputs history field columns as specified by a satellite track data file
           3             : !
           4             : !-------------------------------------------------------------------------------
           5             : module sat_hist
           6             :
           7             :   use perf_mod,            only: t_startf, t_stopf
           8             :   use shr_kind_mod,        only: r4 => shr_kind_r8
           9             :   use shr_kind_mod,        only: r8 => shr_kind_r8, cl=>shr_kind_cl
          10             :   use cam_logfile,         only: iulog
          11             :   use ppgrid,              only: pcols, pver, pverp, begchunk, endchunk
          12             :   use cam_history_support, only: fieldname_lenp2, max_string_len, ptapes
          13             :   use spmd_utils,          only: masterproc, iam
          14             :   use cam_abortutils,      only: endrun
          15             :
          16             :   use pio,                 only: file_desc_t, iosystem_desc_t, var_desc_t, io_desc_t
          17             :   use pio,                 only: pio_inq_dimid, pio_inq_varid
          18             :   use pio,                 only: pio_seterrorhandling, pio_def_var
          19             :   use pio,                 only: pio_inq_dimlen, pio_get_att, pio_put_att, pio_get_var, pio_put_var, pio_write_darray
          20             :   use pio,                 only: pio_real, pio_double
          21             :   use pio,                 only: PIO_NOWRITE, PIO_NOERR, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR, PIO_GLOBAL
          22             :   use spmd_utils,          only: mpicom
          23             : #ifdef SPMD
          24             :   use mpishorthand,        only: mpichar, mpiint
          25             : #endif
          26             :    use physconst,          only: pi
          27             :
          28             :   implicit none
          29             :
          30             :   private
          31             :   save
          32             :
          33             :   public :: sat_hist_readnl
          34             :   public :: sat_hist_init
          35             :   public :: sat_hist_write
          36             :   public :: sat_hist_define
          37             :   public :: is_satfile
          38             :
          39             :   character(len=max_string_len)  :: sathist_track_infile
          40             :   type(file_desc_t) :: infile
          41             :
          42             :   integer :: half_step
          43             :   logical :: has_sat_hist = .false.
          44             :
          45             :   integer :: sathist_nclosest
          46             :   integer :: sathist_ntimestep
          47             :
          48             :   real(r8), allocatable :: obs_lats(:)
          49             :   real(r8), allocatable :: obs_lons(:)
          50             :
          51             :   logical  :: doy_format
          52             :   real(r8) :: first_datetime
          53             :   real(r8) :: last_datetime
          54             :   integer  :: last_start_index
          55             :   integer  :: time_ndx
          56             :   integer  :: t_buffer_size
          57             :   integer, allocatable :: date_buffer(:), time_buffer(:)
          58             :   integer :: sat_tape_num=ptapes-1
          59             :
          60             :
          61             :   ! input file
          62             :   integer :: n_profiles
          63             :   integer :: time_vid, date_vid, lat_vid, lon_vid, instr_vid, orbit_vid, prof_vid, zenith_vid
          64             :
          65             :   integer :: in_julian_vid
          66             :   integer :: in_localtime_vid
          67             :   integer :: in_doy_vid
          68             :   integer :: in_occ_type_vid
          69             :
          70             :   integer :: in_start_col
          71             :
          72             :
          73             :   ! output file
          74             :   type(var_desc_t) :: out_latid, out_lonid, out_dstid, out_instrid, out_zenithid, out_orbid, out_profid
          75             :   type(var_desc_t) :: out_instr_lat_vid, out_instr_lon_vid
          76             :   type(var_desc_t) :: out_obs_date_vid, out_obs_time_vid
          77             :   type(var_desc_t) :: out_julian_vid
          78             :   type(var_desc_t) :: out_localtime_vid
          79             :   type(var_desc_t) :: out_doy_vid
          80             :   type(var_desc_t) :: out_occ_type_vid
          81             :
          82             :   logical, parameter :: debug = .false.
          83             :
          84             :   real(r8), parameter :: rad2deg = 180._r8/pi            ! degrees per radian
          85             :
          86             :   logical :: flds_scanned = .false.
          87             :   logical :: has_phys_srf_flds = .false.
          88             :   logical :: has_phys_lev_flds = .false.
          89             :   logical :: has_phys_ilev_flds = .false.
          90             :   logical :: has_dyn_srf_flds = .false.
          91             :   logical :: has_dyn_lev_flds = .false.
          92             :   logical :: has_dyn_ilev_flds = .false.
          93             :
          94             : contains
          95             :
          96             : !-------------------------------------------------------------------------------
          97             :
          98   129026488 :   logical function is_satfile (file_index)
          99             :     integer, intent(in) :: file_index ! index of file in question
         100   129026488 :     is_satfile = file_index == sat_tape_num
         101   129026488 :   end function is_satfile
         102             :
         103             : !-------------------------------------------------------------------------------
         104        1536 :   subroutine sat_hist_readnl(nlfile, hfilename_spec, mfilt, fincl, nhtfrq, avgflag_pertape)
         105             :
         106             :     use namelist_utils,      only: find_group_name
         107             :     use units,               only: getunit, freeunit
         108             :     use cam_history_support, only: pflds
         109             :     use cam_instance,        only: inst_suffix
         110             :
         111             :     implicit none
         112             :
         113             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         114             :     character(len=*), intent(inout) :: hfilename_spec(:)
         115             :     character(len=*), intent(inout) :: fincl(:,:)
         116             :     character(len=1), intent(inout) :: avgflag_pertape(:)
         117             :     integer,          intent(inout) :: mfilt(:), nhtfrq(:)
         118             :
         119             :     ! Local variables
         120             :     integer :: unitn, ierr
         121             :     character(len=*), parameter :: subname = 'sat_hist_readnl'
         122             :     integer :: f, fcnt
         123             :
         124             :     character(len=fieldname_lenp2) :: sathist_fincl(pflds)
         125             :     character(len=max_string_len)  :: sathist_hfilename_spec
         126             :     integer :: sathist_mfilt, sat_tape_num
         127             :
         128             :     namelist /satellite_options_nl/ sathist_track_infile, sathist_hfilename_spec, sathist_fincl, &
         129             :          sathist_mfilt, sathist_nclosest, sathist_ntimestep
         130             :
         131             :     ! set defaults
         132             :
         133        1536 :     sathist_track_infile = ' '
         134        1536 :     sathist_hfilename_spec = '%c.cam' // trim(inst_suffix) // '.hs.%y-%m-%d-%s.nc'
         135     1537536 :     sathist_fincl(:) = ' '
         136        1536 :     sathist_mfilt = 100000
         137        1536 :     sathist_nclosest = 1
         138        1536 :     sathist_ntimestep = 1
         139             :
         140             :     !read namelist options
         141             :
         142        1536 :     if (masterproc) then
         143           2 :        unitn = getunit()
         144           2 :        open( unitn, file=trim(nlfile), status='old' )
         145           2 :        call find_group_name(unitn, 'satellite_options_nl', status=ierr)
         146           2 :        if (ierr == 0) then
         147           0 :           read(unitn, satellite_options_nl, iostat=ierr)
         148           0 :           if (ierr /= 0) then
         149           0 :              call endrun(subname // ':: ERROR reading namelist')
         150             :           end if
         151             :        end if
         152           2 :        close(unitn)
         153           2 :        call freeunit(unitn)
         154             :     end if
         155             :
         156             : #ifdef SPMD
         157             :     ! broadcast the options to all MPI tasks
         158        1536 :     call mpibcast(sathist_track_infile,   len(sathist_track_infile),   mpichar, 0, mpicom)
         159        1536 :     call mpibcast(sathist_hfilename_spec, len(sathist_hfilename_spec), mpichar, 0, mpicom)
         160        1536 :     call mpibcast(sathist_fincl,          pflds*len(sathist_fincl(1)), mpichar, 0, mpicom)
         161        1536 :     call mpibcast(sathist_mfilt,          1,                           mpiint,  0, mpicom)
         162        1536 :     call mpibcast(sathist_nclosest,       1,                           mpiint,  0, mpicom)
         163        1536 :     call mpibcast(sathist_ntimestep,      1,                           mpiint,  0, mpicom)
         164             : #endif
         165             :
         166        1536 :     has_sat_hist = len_trim(sathist_track_infile) > 0
         167             :
         168        1536 :     if (.not.has_sat_hist) return
         169             :
         170           0 :      sat_tape_num=ptapes-1
         171           0 :      hfilename_spec(sat_tape_num) = sathist_hfilename_spec
         172           0 :      mfilt(sat_tape_num) = sathist_mfilt
         173           0 :      fcnt=0
         174           0 :      do f=1, pflds
         175           0 :         fincl(f,sat_tape_num) = sathist_fincl(f)
         176           0 :         if(len_trim(sathist_fincl(f)) > 0) then
         177           0 :            fcnt=fcnt+1
         178             :         end if
         179             :      enddo
         180             :
         181           0 :      nhtfrq(sat_tape_num) = 1
         182           0 :      avgflag_pertape(sat_tape_num) = 'I'
         183             :
         184           0 :      if(masterproc) then
         185           0 :         write(iulog,*) 'sathist_track_infile: ',trim(sathist_track_infile)
         186           0 :         write(iulog,*) 'sathist_hfilename_spec: ',trim(sathist_hfilename_spec)
         187           0 :         write(iulog,*) 'sathist_fincl: ',(trim(sathist_fincl(f))//' ', f=1,fcnt)
         188           0 :         write(iulog,*) 'max columns per file sathist_mfilt: ',sathist_mfilt
         189           0 :         write(iulog,*) 'sathist_nclosest: ',sathist_nclosest
         190           0 :         write(iulog,*) 'sathist_ntimestep: ',sathist_ntimestep
         191             :      end if
         192             :
         193        1536 :    end subroutine sat_hist_readnl
         194             :
         195             :
         196             : !-------------------------------------------------------------------------------
         197             : !-------------------------------------------------------------------------------
         198         768 :   subroutine sat_hist_init
         199        1536 :     use cam_pio_utils, only: cam_pio_openfile
         200             :     use ioFileMod,     only: getfil
         201             :     use time_manager,  only: get_step_size
         202             :     use string_utils,  only: to_lower, GLC
         203             :
         204             :     implicit none
         205             :
         206             :     character(len=max_string_len)  :: locfn       ! Local filename
         207             :     integer :: ierr, dimid
         208             :
         209             :     character(len=128) :: date_format
         210             :
         211         768 :     if (.not.has_sat_hist) return
         212             :
         213           0 :     call getfil (sathist_track_infile, locfn)
         214           0 :     call cam_pio_openfile(infile, locfn, PIO_NOWRITE)
         215             :
         216           0 :     ierr = pio_inq_dimid(infile,'profs',dimid)
         217           0 :     ierr = pio_inq_dimlen(infile, dimid, n_profiles)
         218             :
         219           0 :     ierr = pio_inq_varid( infile, 'time', time_vid )
         220           0 :     ierr = pio_inq_varid( infile, 'date', date_vid )
         221             :
         222           0 :     ierr = pio_get_att( infile, date_vid, 'long_name', date_format)
         223           0 :     date_format = to_lower(trim( date_format(:GLC(date_format))))
         224             :
         225           0 :     if ( index( date_format, 'yyyymmdd') > 0 ) then
         226           0 :        doy_format = .false.
         227           0 :     else if  ( index( date_format, 'yyyyddd') > 0 ) then
         228           0 :        doy_format = .true.
         229             :     else
         230           0 :        call endrun('sat_hist_init: date_format not recognized : '//trim(date_format))
         231             :     endif
         232             :
         233           0 :     ierr = pio_inq_varid( infile, 'lat', lat_vid )
         234           0 :     ierr = pio_inq_varid( infile, 'lon', lon_vid )
         235             :
         236           0 :     call pio_seterrorhandling(infile, PIO_BCAST_ERROR)
         237           0 :     ierr = pio_inq_varid( infile, 'instr_num', instr_vid )
         238           0 :     if(ierr/=PIO_NOERR) instr_vid=-1
         239             :
         240           0 :     ierr = pio_inq_varid( infile, 'orbit_num', orbit_vid )
         241           0 :     if(ierr/=PIO_NOERR) orbit_vid=-1
         242             :
         243           0 :     ierr = pio_inq_varid( infile, 'prof_num',  prof_vid )
         244           0 :     if(ierr/=PIO_NOERR) prof_vid=-1
         245             :
         246           0 :     ierr = pio_inq_varid( infile, 'instr_sza', zenith_vid )
         247           0 :     if(ierr/=PIO_NOERR) zenith_vid=-1
         248             :
         249           0 :     ierr = pio_inq_varid( infile, 'julian', in_julian_vid )
         250           0 :     if(ierr/=PIO_NOERR) in_julian_vid=-1
         251             :
         252           0 :     ierr = pio_inq_varid( infile, 'local_time', in_localtime_vid )
         253           0 :     if(ierr/=PIO_NOERR) in_localtime_vid=-1
         254             :
         255           0 :     ierr = pio_inq_varid( infile, 'doy', in_doy_vid )
         256           0 :     if(ierr/=PIO_NOERR) in_doy_vid=-1
         257             :
         258           0 :     ierr = pio_inq_varid( infile, 'occ_type', in_occ_type_vid )
         259           0 :     if(ierr/=PIO_NOERR) in_occ_type_vid=-1
         260             :
         261           0 :     call pio_seterrorhandling(infile, PIO_INTERNAL_ERROR)
         262             :
         263           0 :     call read_datetime( first_datetime, 1 )
         264           0 :     call read_datetime( last_datetime, n_profiles )
         265           0 :     last_start_index = -1
         266           0 :     t_buffer_size = min(1000,n_profiles)
         267           0 :     allocate( date_buffer(t_buffer_size), time_buffer(t_buffer_size) )
         268           0 :     if (masterproc) write(iulog,*) "sathist_init:", n_profiles, first_datetime, last_datetime
         269           0 :     if ( last_datetime<first_datetime ) then
         270           0 :        call endrun('sat_hist_init: satellite track file has invalid date time info')
         271             :     endif
         272             :
         273           0 :     time_ndx = 1
         274           0 :     half_step = get_step_size()*0.5_r8
         275             :
         276         768 :   end subroutine sat_hist_init
         277             :
         278             : !-------------------------------------------------------------------------------
         279             : !-------------------------------------------------------------------------------
         280           0 :   subroutine read_datetime( datetime, index )
         281             :
         282             :     real(r8), intent( out ) :: datetime
         283             :     integer,  intent( in )  :: index
         284             :
         285             :     integer :: ierr
         286             :     integer :: cnt(1)
         287             :     integer :: start(1)
         288             :     integer :: date(1), time(1)
         289             :
         290           0 :     cnt = (/ 1 /)
         291           0 :     start = (/index/)
         292             :
         293           0 :     ierr = pio_get_var( infile, time_vid, start, cnt, time )
         294           0 :     ierr = pio_get_var( infile, date_vid, start, cnt, date )
         295             :
         296           0 :     datetime = convert_date_time( date(1),time(1) )
         297             :
         298         768 :   end subroutine read_datetime
         299             :
         300             : !-------------------------------------------------------------------------------
         301             : !-------------------------------------------------------------------------------
         302           0 :   subroutine read_buffered_datetime( datetime, index )
         303             :
         304             :     real(r8), intent( out ) :: datetime
         305             :     integer,  intent( in )  :: index
         306             :
         307             :     integer :: ii
         308             :
         309             :     integer :: ierr
         310             :     integer :: cnt
         311             :     integer :: start
         312             :     integer :: date, time
         313             :
         314             :     ! If the request is outside of the buffer then reload the buffer.
         315             :     if ((last_start_index == -1) .or. (index < last_start_index) &
         316           0 :          .or. (index >= (last_start_index + t_buffer_size))) then
         317             :
         318           0 :        start = (index - 1) / t_buffer_size * t_buffer_size + 1
         319           0 :        if ( start+t_buffer_size-1 <= n_profiles ) then
         320             :           cnt = t_buffer_size
         321             :        else
         322           0 :           cnt = n_profiles-start+1
         323             :        endif
         324           0 :        ierr = pio_get_var( infile, time_vid, (/ start /), (/ cnt /), time_buffer(1:cnt) )
         325           0 :        ierr = pio_get_var( infile, date_vid, (/ start /), (/ cnt /), date_buffer(1:cnt) )
         326             :
         327           0 :        last_start_index = start
         328             :     endif
         329             :
         330           0 :     ii = mod( index - 1, t_buffer_size ) + 1
         331           0 :     time = time_buffer(ii)
         332           0 :     date = date_buffer(ii)
         333           0 :     datetime = convert_date_time( date,time )
         334             :
         335           0 :   end subroutine read_buffered_datetime
         336             :
         337             : !-------------------------------------------------------------------------------
         338             : !-------------------------------------------------------------------------------
         339           0 :   function convert_date_time( date,time )
         340             :     use time_manager, only: set_time_float_from_date
         341             :
         342             :     integer, intent(in) :: date,time
         343             :     real(r8) :: convert_date_time
         344             :
         345             :     real(r8) :: datetime
         346             :     integer :: yr, doy, mon, dom
         347             :
         348           0 :     if ( doy_format ) then
         349           0 :        yr = date/1000
         350           0 :        doy = date - yr*1000
         351           0 :        call set_time_float_from_date( datetime, yr, 1, doy, time )
         352             :     else
         353           0 :        yr = date/10000
         354           0 :        mon = (date - yr*10000)/100
         355           0 :        dom = date - yr*10000 - mon*100
         356           0 :        call set_time_float_from_date( datetime, yr, mon, dom, time )
         357             :     endif
         358           0 :     convert_date_time = datetime
         359             :
         360           0 :   end function convert_date_time
         361             : !-------------------------------------------------------------------------------
         362           0 :   subroutine sat_hist_define(outfile)
         363           0 :     use pio, only : pio_inquire
         364             :     type(file_desc_t), intent(inout) :: outfile
         365             :
         366             :     integer :: coldim
         367             :     integer :: ierr
         368             :
         369           0 :     ierr = pio_inquire(outfile, unlimitedDimId=coldim)
         370             :
         371           0 :     call pio_seterrorhandling(outfile, PIO_BCAST_ERROR)
         372           0 :     ierr = define_var( 'instr_lat', coldim, infile, lat_vid,  outfile, out_instr_lat_vid )
         373           0 :     ierr = define_var( 'instr_lon', coldim, infile, lon_vid,  outfile, out_instr_lon_vid )
         374           0 :     ierr = define_var( 'obs_time', coldim, infile, time_vid,  outfile, out_obs_time_vid )
         375           0 :     ierr = define_var( 'obs_date', coldim, infile, date_vid,  outfile, out_obs_date_vid )
         376             :
         377           0 :     ierr = pio_inq_varid( outfile, 'distance', out_dstid )
         378           0 :     if (ierr /= PIO_NOERR) then
         379           0 :        ierr = pio_def_var  ( outfile, 'distance', PIO_REAL, (/coldim/), out_dstid )
         380           0 :        ierr = pio_put_att  ( outfile, out_dstid, "long_name", "distance from midpoint to observation")
         381           0 :        ierr = pio_put_att  ( outfile, out_dstid, "units", "km")
         382             :     end if
         383             :
         384           0 :     if (orbit_vid>0) then
         385           0 :        ierr = define_var( 'orbit_num', coldim, infile, orbit_vid,  outfile, out_orbid )
         386             :     endif
         387           0 :     if (prof_vid>0) then
         388           0 :        ierr = define_var( 'prof_num', coldim, infile, prof_vid,  outfile, out_profid )
         389             :     endif
         390           0 :     if (instr_vid>0) then
         391           0 :        ierr = define_var( 'instr_num', coldim, infile, instr_vid,  outfile, out_instrid )
         392             :     endif
         393           0 :     if (zenith_vid>0) then
         394           0 :        ierr = define_var( 'instr_sza', coldim, infile, zenith_vid,  outfile, out_zenithid )
         395             :     endif
         396           0 :     if (in_occ_type_vid>0) then
         397           0 :        ierr = define_var( 'occ_type', coldim, infile, in_occ_type_vid,  outfile, out_occ_type_vid )
         398             :     endif
         399           0 :     if (in_julian_vid>0) then
         400           0 :        ierr = define_var( 'julian', coldim, infile, in_julian_vid,  outfile, out_julian_vid )
         401             :     endif
         402           0 :     if (in_localtime_vid>0) then
         403           0 :        ierr = define_var( 'local_time', coldim, infile, in_localtime_vid,  outfile, out_localtime_vid )
         404             :     endif
         405           0 :     if (in_doy_vid>0) then
         406           0 :        ierr = define_var( 'doy', coldim, infile, in_doy_vid,  outfile, out_doy_vid )
         407             :     endif
         408             :
         409           0 :     call pio_seterrorhandling(outfile, PIO_INTERNAL_ERROR)
         410           0 :     ierr=pio_put_att (outfile, PIO_GLOBAL, 'satellite_track_file', sathist_track_infile)
         411           0 :   end subroutine sat_hist_define
         412             :
         413             : !-------------------------------------------------------------------------------
         414           0 :   subroutine sat_hist_write( tape , nflds, nfils)
         415             :
         416             :     use phys_grid, only: phys_decomp
         417             :     use dyn_grid,  only: dyn_decomp
         418             :     use cam_history_support, only : active_entry
         419             :     use pio, only : pio_file_is_open, pio_syncfile
         420             :
         421             :     type(active_entry) :: tape
         422             :     integer, intent(in) :: nflds
         423             :     integer, intent(inout) :: nfils
         424             :
         425             :     integer :: ncols, nocols
         426             :     integer :: ierr
         427             :
         428           0 :     integer, allocatable :: col_ndxs(:)
         429           0 :     integer, allocatable :: chk_ndxs(:)
         430           0 :     integer, allocatable :: fdyn_ndxs(:)
         431           0 :     integer, allocatable :: ldyn_ndxs(:)
         432           0 :     integer, allocatable :: phs_owners(:)
         433           0 :     integer, allocatable :: dyn_owners(:)
         434           0 :     real(r8),allocatable :: mlats(:)
         435           0 :     real(r8),allocatable :: mlons(:)
         436           0 :     real(r8),allocatable :: phs_dists(:)
         437             :
         438             :     integer :: coldim
         439             :     logical :: has_dyn_flds
         440             :
         441           0 :     if (.not.has_sat_hist) return
         442             :
         443           0 :     call read_next_position( ncols )
         444             :
         445           0 :     if ( ncols < 1 ) return
         446             :
         447           0 :     call t_startf ('sat_hist_write')
         448             :
         449             :     ! The n closest columns to the observation will be output,
         450             :     ! so increase the size of the columns used for output/
         451           0 :     nocols = ncols * sathist_nclosest
         452             :
         453           0 :     allocate( col_ndxs(nocols) )
         454           0 :     allocate( chk_ndxs(nocols) )
         455           0 :     allocate( fdyn_ndxs(nocols) )
         456           0 :     allocate( ldyn_ndxs(nocols) )
         457           0 :     allocate( phs_owners(nocols) )
         458           0 :     allocate( dyn_owners(nocols) )
         459           0 :     allocate( mlats(nocols) )
         460           0 :     allocate( mlons(nocols) )
         461           0 :     allocate( phs_dists(nocols) )
         462             :
         463           0 :     call scan_flds( tape, nflds )
         464           0 :     has_dyn_flds = has_dyn_srf_flds .or. has_dyn_lev_flds .or. has_dyn_ilev_flds
         465             :
         466             :     call get_indices( obs_lats, obs_lons, ncols, nocols, has_dyn_flds, col_ndxs, chk_ndxs, &
         467           0 :          fdyn_ndxs, ldyn_ndxs, phs_owners, dyn_owners, mlats, mlons, phs_dists )
         468             :
         469           0 :     if ( .not. pio_file_is_open(tape%Files(1)) ) then
         470           0 :        call endrun('sat file not open')
         471             :     endif
         472             :
         473             :
         474           0 :     ierr = pio_inq_dimid(tape%Files(1),'ncol',coldim )
         475             :
         476           0 :     ierr = pio_inq_varid(tape%Files(1), 'lat', out_latid )
         477           0 :     ierr = pio_inq_varid(tape%Files(1), 'lon', out_lonid )
         478           0 :     ierr = pio_inq_varid(tape%Files(1), 'distance', out_dstid )
         479             :
         480           0 :     call write_record_coord( tape, mlats(:), mlons(:), phs_dists(:), ncols, nfils )
         481             :
         482             :    ! dump columns of 2D fields
         483           0 :     if (has_phys_srf_flds) then
         484             :        call dump_columns( tape%Files(1), tape%hlist, nflds, nocols, 1, nfils, &
         485           0 :                           col_ndxs, chk_ndxs, phs_owners, phys_decomp )
         486             :     endif
         487           0 :     if (has_dyn_srf_flds) then
         488             :        call dump_columns( tape%Files(1), tape%hlist, nflds, nocols, 1, nfils, &
         489           0 :                           fdyn_ndxs, ldyn_ndxs, dyn_owners, dyn_decomp )
         490             :     endif
         491             :
         492             :    ! dump columns of 3D fields defined on mid pres levels
         493           0 :     if (has_phys_lev_flds) then
         494             :        call dump_columns( tape%Files(1), tape%hlist, nflds, nocols, pver, nfils, &
         495           0 :                           col_ndxs, chk_ndxs, phs_owners, phys_decomp )
         496             :     endif
         497           0 :     if (has_dyn_lev_flds) then
         498             :        call dump_columns( tape%Files(1), tape%hlist, nflds, nocols, pver, nfils, &
         499           0 :                           fdyn_ndxs, ldyn_ndxs, dyn_owners, dyn_decomp )
         500             :     endif
         501             :
         502             :    ! dump columns of 3D fields defined on interface pres levels
         503           0 :     if (has_phys_ilev_flds) then
         504             :        call dump_columns( tape%Files(1), tape%hlist, nflds, nocols, pverp, nfils, &
         505           0 :                           col_ndxs, chk_ndxs, phs_owners, phys_decomp )
         506             :     endif
         507           0 :     if (has_dyn_ilev_flds) then
         508             :        call dump_columns( tape%Files(1), tape%hlist, nflds, nocols, pverp, nfils, &
         509           0 :                           fdyn_ndxs, ldyn_ndxs, dyn_owners, dyn_decomp )
         510             :     endif
         511             :
         512           0 :     deallocate( col_ndxs, chk_ndxs, fdyn_ndxs, ldyn_ndxs, phs_owners, dyn_owners )
         513           0 :     deallocate( mlons, mlats, phs_dists )
         514           0 :     deallocate( obs_lons, obs_lats )
         515           0 :     call pio_syncfile(tape%Files(1))
         516             :
         517           0 :     nfils = nfils + nocols
         518             :
         519           0 :     call t_stopf ('sat_hist_write')
         520             :
         521           0 :   end subroutine sat_hist_write
         522             :
         523             : !-------------------------------------------------------------------------------
         524           0 :   subroutine dump_columns( File, hitems, nflds, ncols, nlevs, nfils, fdims, ldims, owners, decomp )
         525           0 :     use cam_history_support, only: field_info, hentry, fillvalue
         526             :     use pio,                 only: pio_setframe, pio_offset_kind
         527             :     use spmd_utils, only:  mpi_real4, mpi_real8, mpicom, mpi_sum
         528             :
         529             :     type(File_desc_t),intent(inout)  :: File
         530             :     type(hentry),     intent(in), target :: hitems(:)
         531             :     integer,          intent(in)     :: nflds
         532             :     integer,          intent(in)     :: ncols
         533             :     integer,          intent(in)     :: nlevs
         534             :     integer,          intent(in)     :: nfils
         535             :     integer,          intent(in)     :: fdims(:)
         536             :     integer,          intent(in)     :: ldims(:)
         537             :     integer,          intent(in)     :: owners(:)
         538             :     integer,          intent(in)     :: decomp
         539             :
         540             :
         541             :     type(field_info), pointer :: field
         542             :     type(var_desc_t) :: vardesc
         543             :     integer :: ierr
         544             :
         545           0 :     real(r8) :: sbuf1d(ncols),rbuf1d(ncols)
         546           0 :     real(r4) :: buf1d(ncols)
         547           0 :     real(r8) :: sbuf2d(nlevs,ncols), rbuf2d(nlevs,ncols)
         548           0 :     real(r4) :: buf2d(nlevs,ncols)
         549             :     integer :: i,k,f, cnt
         550             :
         551           0 :     call t_startf ('sat_hist::dump_columns')
         552             :
         553           0 :     do f = 1,nflds
         554           0 :        field => hitems(f)%field
         555             :
         556           0 :        if (field%numlev==nlevs .and. field%decomp_type==decomp) then
         557           0 :           vardesc = hitems(f)%varid(1)
         558             :
         559           0 :           if (nlevs==1) then
         560           0 :              sbuf1d = 0.0_r8
         561           0 :              rbuf1d = 0.0_r8
         562           0 :              do i=1,ncols
         563           0 :                 if ( iam == owners(i) ) then
         564           0 :                    sbuf1d(i) = hitems(f)%hbuf( fdims(i), 1, ldims(i) )
         565             :                 endif
         566             :              enddo
         567           0 :              call mpi_allreduce(sbuf1d,rbuf1d,ncols,mpi_real8, mpi_sum, mpicom, ierr)
         568           0 :              buf1d(:) = real(rbuf1d(:),r4)
         569           0 :              ierr = pio_put_var(File, vardesc, (/nfils/),(/ncols/), buf1d(:))
         570             :           else
         571           0 :              sbuf2d = 0.0_r8
         572           0 :              rbuf2d = 0.0_r8
         573           0 :              do i=1,ncols
         574           0 :                 if ( iam == owners(i) ) then
         575           0 :                    do k = 1,nlevs
         576           0 :                       sbuf2d(k,i) = hitems(f)%hbuf( fdims(i), k, ldims(i) )
         577             :                    enddo
         578             :                 endif
         579             :              enddo
         580           0 :              call mpi_allreduce(sbuf2d,rbuf2d,ncols*nlevs,mpi_real8, mpi_sum, mpicom, ierr)
         581           0 :              buf2d(:,:) = real(rbuf2d(:,:),r4)
         582           0 :              ierr = pio_put_var(File, vardesc, (/1,nfils/),(/nlevs,ncols/), buf2d(:,:))
         583             :           endif
         584             :
         585             :        endif
         586             :
         587             :     enddo
         588             :
         589           0 :     call t_stopf ('sat_hist::dump_columns')
         590             :
         591           0 :   end subroutine dump_columns
         592             :
         593             : !-------------------------------------------------------------------------------
         594             : ! scan the fields for possible different decompositions
         595             : !-------------------------------------------------------------------------------
         596           0 :   subroutine scan_flds( tape, nflds )
         597           0 :     use cam_history_support, only : active_entry
         598             :     use phys_grid, only: phys_decomp
         599             :     use dyn_grid,  only: dyn_decomp
         600             :
         601             :     type(active_entry), intent(in) :: tape
         602             :     integer, intent(in) :: nflds
         603             :
         604             :     integer :: f
         605             :     character(len=cl) :: msg1, msg2
         606             :
         607           0 :     if (flds_scanned) return
         608             :
         609           0 :     do f = 1,nflds
         610           0 :        if ( tape%hlist(f)%field%decomp_type == phys_decomp ) then
         611           0 :           if ( tape%hlist(f)%field%numlev == 1 ) then
         612           0 :              has_phys_srf_flds = .true.
         613           0 :           elseif ( tape%hlist(f)%field%numlev == pver ) then
         614           0 :              has_phys_lev_flds = .true.
         615           0 :           elseif ( tape%hlist(f)%field%numlev == pverp ) then
         616           0 :              has_phys_ilev_flds = .true.
         617             :           else
         618           0 :              call endrun('sat_hist::scan_flds numlev error : '//tape%hlist(f)%field%name)
         619             :           endif
         620           0 :        elseif ( tape%hlist(f)%field%decomp_type == dyn_decomp ) then
         621           0 :           if ( tape%hlist(f)%field%numlev == 1 ) then
         622           0 :              has_dyn_srf_flds = .true.
         623           0 :           elseif ( tape%hlist(f)%field%numlev == pver ) then
         624           0 :              has_dyn_lev_flds = .true.
         625           0 :           elseif ( tape%hlist(f)%field%numlev == pverp ) then
         626           0 :              has_dyn_ilev_flds = .true.
         627             :           else
         628           0 :              call endrun('sat_hist::scan_flds numlev error : '//tape%hlist(f)%field%name)
         629             :           endif
         630             :        else
         631           0 :           call endrun('sat_hist::scan_flds decomp_type error : '//tape%hlist(f)%field%name)
         632             :        endif
         633             :
         634             :        ! Check that the only "mdim" is the vertical coordinate.
         635             :        if (has_phys_srf_flds .or. has_phys_lev_flds .or. has_phys_ilev_flds .or. &
         636           0 :            has_dyn_srf_flds .or. has_dyn_lev_flds .or. has_dyn_ilev_flds) then
         637             :           ! The mdims pointer is unassociated on a restart.  The restart initialization
         638             :           ! should be fixed rather than requiring the check to make sure it is associated.
         639           0 :           if (associated(tape%hlist(f)%field%mdims)) then
         640           0 :              if (size(tape%hlist(f)%field%mdims) > 1) then
         641           0 :                 msg1 = 'sat_hist::scan_flds mdims error :'//tape%hlist(f)%field%name
         642             :                 msg2 = trim(msg1)//' has mdims in addition to the vertical coordinate.'//&
         643           0 :                    new_line('a')//'  This is not currently supported.'
         644           0 :                 write(iulog,*) msg2
         645           0 :                 call endrun(msg1)
         646             :              end if
         647             :           end if
         648             :        end if
         649             :
         650             :     enddo
         651             :
         652           0 :     flds_scanned = .true.
         653           0 :   end subroutine scan_flds
         654             :
         655             : !-------------------------------------------------------------------------------
         656             : !-------------------------------------------------------------------------------
         657           0 :   subroutine read_next_position( ncols )
         658           0 :     use time_manager, only: get_curr_date
         659             :     use time_manager, only: set_time_float_from_date
         660             :
         661             :     implicit none
         662             :
         663             :     integer,  intent(out) :: ncols
         664             :
         665             :     integer :: ierr
         666             :     integer :: yr, mon, day, tod
         667             :     real(r8) :: begdatetime, enddatetime
         668             :     integer :: beg_ndx, end_ndx, i
         669             :
         670             :     real(r8) :: datetime
         671             :
         672           0 :     call get_curr_date(yr, mon, day, tod)
         673           0 :     call set_time_float_from_date(begdatetime, yr, mon, day, tod-half_step*sathist_ntimestep)
         674           0 :     call set_time_float_from_date(enddatetime, yr, mon, day, tod+half_step*sathist_ntimestep)
         675             :
         676           0 :     ncols = 0
         677             :
         678           0 :     if ( first_datetime > enddatetime ) then
         679           0 :        if (masterproc) write(iulog,'(a,2f16.6)') &
         680           0 :             'sat_hist->read_next_position: all of the satellite date times are after the time window', first_datetime, enddatetime
         681           0 :        return
         682             :     endif
         683           0 :     if ( last_datetime < begdatetime ) then
         684           0 :        if (masterproc) write(iulog,'(a,2f16.6)') &
         685           0 :             'sat_hist->read_next_position: all of the satellite date times are before the time window', begdatetime, last_datetime
         686           0 :        return
         687             :     endif
         688             :
         689           0 :     call t_startf ('sat_hist::read_next_position')
         690             :
         691           0 :     beg_ndx = -99
         692           0 :     end_ndx = -99
         693             :
         694           0 :     bnds_loop: do i = time_ndx,n_profiles
         695             :
         696           0 :        call read_buffered_datetime( datetime, i )
         697             :
         698           0 :        if ( datetime>begdatetime .and. beg_ndx<0 ) beg_ndx = i
         699           0 :        if ( datetime>enddatetime ) exit bnds_loop
         700           0 :        end_ndx = i
         701             :     enddo bnds_loop
         702             :
         703           0 :     if (beg_ndx == -99 .and. end_ndx== -99) then
         704           0 :        if (masterproc) write(iulog,'(a)')  'sat_hist->read_next_position: must be beyond last position -- returning.'
         705           0 :        return
         706             :     endif
         707             :
         708             :     ! Advance the search forward, but because of ntimesteps, it is possible
         709             :     ! for observations used here to be used again. However, we should not go
         710             :     ! back before the previous beginning time.
         711           0 :     if (beg_ndx>0) time_ndx = beg_ndx
         712             :
         713           0 :     ncols = end_ndx-beg_ndx+1
         714           0 :     if (ncols > 0) then
         715           0 :        allocate( obs_lats(ncols), obs_lons(ncols) )
         716           0 :        in_start_col = beg_ndx
         717             :
         718           0 :        ierr = pio_get_var( infile, lat_vid, (/beg_ndx/), (/ncols/), obs_lats )
         719           0 :        ierr = pio_get_var( infile, lon_vid, (/beg_ndx/), (/ncols/), obs_lons )
         720             :
         721             :     endif
         722             :
         723           0 :     call t_stopf ('sat_hist::read_next_position')
         724           0 :   end subroutine read_next_position
         725             :
         726             : !-------------------------------------------------------------------------------
         727             : !-------------------------------------------------------------------------------
         728           0 :   subroutine write_record_coord( tape, mod_lats, mod_lons, mod_dists, ncols, nfils )
         729             :
         730           0 :     use time_manager, only: get_curr_date, get_curr_time
         731             :     use cam_history_support, only : active_entry
         732             :     implicit none
         733             :     type(active_entry), intent(inout) :: tape
         734             :
         735             :     integer,  intent(in) :: ncols
         736             :     real(r8), intent(in) :: mod_lats(ncols * sathist_nclosest)
         737             :     real(r8), intent(in) :: mod_lons(ncols * sathist_nclosest)
         738             :     real(r8), intent(in) :: mod_dists(ncols * sathist_nclosest)
         739             :     integer,  intent(in) :: nfils
         740             :
         741             :     integer :: ierr, i
         742             :     integer :: yr, mon, day      ! year, month, and day components of a date
         743             :     integer :: ncdate            ! current date in integer format [yyyymmdd]
         744             :     integer :: ncsec             ! current time of day [seconds]
         745             :     integer :: ndcur             ! day component of current time
         746             :     integer :: nscur             ! seconds component of current time
         747             :     real(r8) :: time             ! current time
         748           0 :     integer, allocatable  :: itmp(:)
         749           0 :     real(r8), allocatable :: rtmp(:)
         750           0 :     real(r8), allocatable :: out_lats(:)
         751           0 :     real(r8), allocatable :: out_lons(:)
         752             :
         753           0 :     call t_startf ('sat_hist::write_record_coord')
         754             :
         755           0 :     call get_curr_date(yr, mon, day, ncsec)
         756           0 :     ncdate = yr*10000 + mon*100 + day
         757           0 :     call get_curr_time(ndcur, nscur)
         758             :
         759             :
         760           0 :     time = ndcur + nscur/86400._r8
         761             :
         762           0 :     allocate( itmp(ncols * sathist_nclosest) )
         763           0 :     allocate( rtmp(ncols * sathist_nclosest) )
         764             :
         765           0 :     itmp(:) = ncdate
         766           0 :     ierr = pio_put_var(tape%Files(1), tape%dateid,(/nfils/), (/ncols * sathist_nclosest/),itmp)
         767           0 :     itmp(:) = ncsec
         768           0 :     ierr = pio_put_var(tape%Files(1), tape%datesecid,(/nfils/),(/ncols * sathist_nclosest/),itmp)
         769           0 :     rtmp(:) = time
         770           0 :     ierr = pio_put_var(tape%Files(1), tape%timeid, (/nfils/),(/ncols * sathist_nclosest/),rtmp)
         771             :
         772           0 :     deallocate(itmp)
         773           0 :     deallocate(rtmp)
         774             :
         775             :     ! output model column coordinates
         776           0 :     ierr = pio_put_var(tape%Files(1), out_latid, (/nfils/),(/ncols * sathist_nclosest/), mod_lats)
         777           0 :     ierr = pio_put_var(tape%Files(1), out_lonid, (/nfils/),(/ncols * sathist_nclosest/), mod_lons)
         778           0 :     ierr = pio_put_var(tape%Files(1), out_dstid, (/nfils/),(/ncols * sathist_nclosest/), mod_dists / 1000._r8)
         779             :
         780             :     ! output instrument location
         781           0 :     allocate( out_lats(ncols * sathist_nclosest) )
         782           0 :     allocate( out_lons(ncols * sathist_nclosest) )
         783             :
         784           0 :     do i = 1, ncols
         785           0 :       out_lats(((i-1)*sathist_nclosest)+1 : (i*sathist_nclosest)) = obs_lats(i)
         786           0 :       out_lons(((i-1)*sathist_nclosest)+1 : (i*sathist_nclosest)) = obs_lons(i)
         787             :     enddo
         788             :
         789           0 :     ierr = pio_put_var(tape%Files(1), out_instr_lat_vid, (/nfils/),(/ncols * sathist_nclosest/), out_lats)
         790           0 :     ierr = pio_put_var(tape%Files(1), out_instr_lon_vid, (/nfils/),(/ncols * sathist_nclosest/), out_lons)
         791             :
         792           0 :     deallocate(out_lats)
         793           0 :     deallocate(out_lons)
         794             :
         795             :
         796           0 :     ierr = copy_data( infile, date_vid, tape%Files(1), out_obs_date_vid, in_start_col, nfils, ncols )
         797           0 :     ierr = copy_data( infile, time_vid, tape%Files(1), out_obs_time_vid, in_start_col, nfils, ncols )
         798             :
         799             :     ! output observation identifiers
         800           0 :     if (instr_vid>0) then
         801           0 :        ierr = copy_data( infile, instr_vid, tape%Files(1), out_instrid, in_start_col, nfils, ncols )
         802             :     endif
         803           0 :     if (orbit_vid>0) then
         804           0 :        ierr = copy_data( infile, orbit_vid, tape%Files(1), out_orbid, in_start_col, nfils, ncols )
         805             :     endif
         806           0 :     if (prof_vid>0) then
         807           0 :        ierr = copy_data( infile, prof_vid, tape%Files(1), out_profid, in_start_col, nfils, ncols )
         808             :     endif
         809           0 :     if (zenith_vid>0) then
         810           0 :        ierr = copy_data( infile, zenith_vid, tape%Files(1), out_zenithid, in_start_col, nfils, ncols )
         811             :     endif
         812           0 :     if (in_julian_vid>0) then
         813           0 :        ierr = copy_data( infile, in_julian_vid, tape%Files(1), out_julian_vid, in_start_col, nfils, ncols )
         814             :     endif
         815           0 :     if (in_occ_type_vid>0) then
         816           0 :        ierr = copy_data( infile, in_occ_type_vid, tape%Files(1), out_occ_type_vid, in_start_col, nfils, ncols )
         817             :     endif
         818           0 :     if (in_localtime_vid>0) then
         819           0 :        ierr = copy_data( infile, in_localtime_vid, tape%Files(1), out_localtime_vid, in_start_col, nfils, ncols )
         820             :     endif
         821           0 :     if (in_doy_vid>0) then
         822           0 :        ierr = copy_data( infile, in_doy_vid, tape%Files(1), out_doy_vid, in_start_col, nfils, ncols )
         823             :     endif
         824             :
         825           0 :     call t_stopf ('sat_hist::write_record_coord')
         826           0 :   end subroutine write_record_coord
         827             :
         828             : !-------------------------------------------------------------------------------
         829             : !-------------------------------------------------------------------------------
         830             :
         831           0 :   subroutine get_indices( lats, lons, ncols, nocols, has_dyn_flds, col_ndxs, chk_ndxs, &
         832           0 :        fdyn_ndxs, ldyn_ndxs, phs_owners, dyn_owners, mlats, mlons, phs_dists )
         833             :
         834           0 :     use dyn_grid, only : dyn_grid_get_colndx
         835             :
         836             :     integer,  intent(in)  :: ncols
         837             :     real(r8), intent(in)  :: lats(ncols)
         838             :     real(r8), intent(in)  :: lons(ncols)
         839             :     integer,  intent(in)  :: nocols
         840             :     logical,  intent(in)  :: has_dyn_flds
         841             :     integer,  intent(out) :: col_ndxs(nocols)
         842             :     integer,  intent(out) :: chk_ndxs(nocols)
         843             :     integer,  intent(out) :: fdyn_ndxs(nocols)
         844             :     integer,  intent(out) :: ldyn_ndxs(nocols)
         845             :     integer,  intent(out) :: phs_owners(nocols)
         846             :     integer,  intent(out) :: dyn_owners(nocols)
         847             :     real(r8), intent(out) :: mlats(nocols)
         848             :     real(r8), intent(out) :: mlons(nocols)
         849             :     real(r8), intent(out) :: phs_dists(nocols)
         850             :
         851             :     integer :: i, j, ndx
         852             :     real(r8) :: lat, lon
         853             :
         854           0 :     integer,  allocatable :: ichks(:),icols(:),idyn1s(:),idyn2s(:), iphs_owners(:), idyn_owners(:)
         855           0 :     real(r8), allocatable :: rlats(:), rlons(:), plats(:), plons(:), iphs_dists(:)
         856             :
         857           0 :     integer :: gcols(sathist_nclosest)
         858             :
         859           0 :     call t_startf ('sat_hist::get_indices')
         860             :
         861             :     allocate(ichks(sathist_nclosest),icols(sathist_nclosest),idyn1s(sathist_nclosest), &
         862           0 :          idyn2s(sathist_nclosest),iphs_owners(sathist_nclosest),idyn_owners(sathist_nclosest))
         863             :     allocate(rlats(sathist_nclosest), rlons(sathist_nclosest), plats(sathist_nclosest), &
         864           0 :          plons(sathist_nclosest), iphs_dists(sathist_nclosest) )
         865             :
         866           0 :     col_ndxs = -1
         867           0 :     chk_ndxs = -1
         868           0 :     fdyn_ndxs = -1
         869           0 :     ldyn_ndxs = -1
         870           0 :     phs_owners = -1
         871           0 :     dyn_owners = -1
         872           0 :     phs_dists = -1
         873             :
         874             :     ndx = 0
         875           0 :     do i = 1,ncols
         876             :
         877           0 :        lat = lats(i)
         878           0 :        lon = lons(i)
         879             :
         880           0 :        if ( lon >= 360._r8) then
         881           0 :          lon = lon-360._r8
         882             :        endif
         883           0 :        if ( lon < 0._r8) then
         884           0 :          lon = lon+360._r8
         885             :        endif
         886           0 :        if (lat<-90._r8 .or. lat>90._r8) then
         887           0 :           write(iulog,*) 'sat_hist::get_indices lat = ',lat
         888           0 :           call endrun('sat_hist::get_indices : lat must be between -90 and 90 degrees (-90<=lat<=90)')
         889             :        endif
         890           0 :        if (lon<0._r8 .or. lon>=360._r8) then
         891           0 :           write(iulog,*) 'sat_hist::get_indices lon = ',lon
         892           0 :           call endrun('sat_hist::get_indices : lon must be between 0 and 360 degrees (0<=lon<360)')
         893             :        endif
         894             :
         895             :        call find_cols( lat, lon, sathist_nclosest, iphs_owners, ichks, icols, &
         896           0 :                        gcols, iphs_dists, plats, plons )
         897             :
         898           0 :        if (has_dyn_flds) then
         899           0 :           call dyn_grid_get_colndx( gcols, sathist_nclosest, idyn_owners, idyn1s, idyn2s )
         900             :        endif
         901             :
         902           0 :        do j = 1, sathist_nclosest
         903             :
         904           0 :           if (debug .and. iam==iphs_owners(j) ) then
         905             :              if ( abs(plats(j)-rlats(j))>1.e-3_r8 ) then
         906             :                 write(*,'(a,3f20.12)') ' lat, plat, rlat = ', lat, plats(j), rlats(j)
         907             :                 write(*,'(a,3f20.12)') ' lon, plon, rlon = ', lon, plons(j), rlons(j)
         908             :                 call endrun('sat_hist::get_indices: dyn lat is different than phys lat ')
         909             :              endif
         910             :              if ( abs(plons(j)-rlons(j))>1.e-3_r8 ) then
         911             :                 write(*,'(a,3f20.12)') ' lat, plat, rlat = ', lat, plats(j), rlats(j)
         912             :                 write(*,'(a,3f20.12)') ' lon, plon, rlon = ', lon, plons(j), rlons(j)
         913             :                 call endrun('sat_hist::get_indices: dyn lon is different than phys lon ')
         914             :              endif
         915             :           endif
         916             :
         917           0 :           ndx = ndx+1
         918             :
         919           0 :           chk_ndxs(ndx)   = ichks(j)
         920           0 :           col_ndxs(ndx)   = icols(j)
         921           0 :           fdyn_ndxs(ndx)  = idyn1s(j)
         922           0 :           ldyn_ndxs(ndx)  = idyn2s(j)
         923           0 :           mlats(ndx)      = plats(j)
         924           0 :           mlons(ndx)      = plons(j)
         925           0 :           phs_owners(ndx) = iphs_owners(j)
         926           0 :           dyn_owners(ndx) = idyn_owners(j)
         927           0 :           phs_dists(ndx)  = iphs_dists(j)
         928             :        enddo
         929             :     enddo
         930             :
         931           0 :     deallocate(ichks, icols, idyn1s, idyn2s, iphs_owners, idyn_owners)
         932           0 :     deallocate(rlats, rlons, plats, plons, iphs_dists )
         933             :
         934           0 :     call t_stopf ('sat_hist::get_indices')
         935           0 :   end subroutine get_indices
         936             :
         937             : !-------------------------------------------------------------------------------
         938             : ! utility function
         939             : !-------------------------------------------------------------------------------
         940           0 :   integer function define_var( var_name, coldim, infile, in_vid, outfile, out_id ) result(res)
         941             :
         942           0 :     use pio, only: pio_inq_vartype
         943             :
         944             :     character(len=*), intent(in) :: var_name
         945             :     integer,          intent(in) :: coldim
         946             :     type(File_desc_t),intent(inout) :: infile
         947             :     type(File_desc_t),intent(inout) :: outfile
         948             :     integer,          intent(in) :: in_vid
         949             :     type(var_desc_t), intent(out):: out_id
         950             :
         951             :     integer :: type
         952             :
         953           0 :     res = pio_inq_varid( outfile, var_name, out_id )
         954           0 :     if(res/=PIO_NOERR) then
         955             :
         956           0 :        res = pio_inq_vartype( infile, in_vid, type )
         957             :
         958           0 :        res = pio_def_var ( outfile, var_name, type, (/coldim/), out_id )
         959             :
         960           0 :        res = copy_att( infile, in_vid, 'long_name', outfile, out_id )
         961           0 :        res = copy_att( infile, in_vid, 'units',     outfile, out_id )
         962             :
         963             :     endif
         964             :
         965           0 :   end function define_var
         966             :
         967             : !-------------------------------------------------------------------------------
         968             : ! utility function
         969             : !-------------------------------------------------------------------------------
         970           0 :   integer function copy_data( infile, in_vid, outfile, out_id, instart, outstart, ncols ) result(res)
         971             :
         972             :     type(File_desc_t),intent(in) :: infile
         973             :     type(File_desc_t),intent(inout) :: outfile
         974             :     integer,          intent(in) :: in_vid
         975             :     type(var_desc_t), intent(in) :: out_id
         976             :     integer,          intent(in) :: instart, outstart, ncols
         977             :
         978           0 :     real(r8), allocatable :: data(:)
         979           0 :     real(r8), allocatable :: outdata(:)
         980             :     integer               :: i
         981             :
         982           0 :     allocate( data(ncols) )
         983             :
         984           0 :     res = pio_get_var( infile,  in_vid, (/instart/),  (/ncols/), data )
         985             :
         986           0 :     allocate( outdata(ncols * sathist_nclosest) )
         987             :
         988           0 :     do i = 1, ncols
         989           0 :       outdata(((i-1)*sathist_nclosest)+1 : (i*sathist_nclosest)) = data(i)
         990             :     enddo
         991             :
         992           0 :     res = pio_put_var( outfile, out_id, (/outstart/), (/ncols * sathist_nclosest/), outdata )
         993             :
         994           0 :     deallocate(outdata)
         995           0 :     deallocate(data)
         996             :
         997           0 :   end function copy_data
         998             :
         999             : !-------------------------------------------------------------------------------
        1000             : ! utility function
        1001             : ! -- should be able to use pio_copy_att which does not seem to work
        1002             : !-------------------------------------------------------------------------------
        1003           0 :   integer function copy_att( infile, in_vid, att_name, outfile, out_id ) result(res)
        1004             :
        1005             :     type(File_desc_t),intent(inout) :: infile
        1006             :     type(File_desc_t),intent(inout) :: outfile
        1007             :     character(len=*), intent(in) :: att_name
        1008             :     integer,          intent(in) :: in_vid
        1009             :     type(var_desc_t), intent(in) :: out_id
        1010             :
        1011             :     character(len=1024) :: att
        1012             :
        1013             :
        1014           0 :     res = pio_get_att( infile, in_vid, trim(att_name), att )
        1015           0 :     if (res==PIO_NOERR) then
        1016           0 :        res = pio_put_att ( outfile, out_id, trim(att_name), trim(att))
        1017             :     endif
        1018             :
        1019             :
        1020           0 :   end function copy_att
        1021             :
        1022             :   !-------------------------------------------------------------------------------
        1023             :   !-------------------------------------------------------------------------------
        1024           0 :   subroutine find_cols(lat, lon, nclosest, owner, lcid, icol, gcol, distmin, mlats, mlons)
        1025             :     use physconst,  only: rearth
        1026             :     use phys_grid,  only: get_rlon_all_p, get_rlat_all_p, get_gcol_p, get_ncols_p
        1027             :     use spmd_utils, only: iam, npes, mpi_integer, mpi_real8, mpicom
        1028             :
        1029             :     real(r8),intent(in)  :: lat, lon            ! requested location in degrees
        1030             :     integer, intent(in)  :: nclosest            ! number of closest points to find
        1031             :     integer, intent(out) :: owner(nclosest)     ! rank of chunk owner
        1032             :     integer, intent(out) :: lcid(nclosest)      ! local chunk index
        1033             :     integer, intent(out) :: icol(nclosest)      ! column index within the chunk
        1034             :     integer, intent(out) :: gcol(nclosest)      ! global column index
        1035             :     real(r8),intent(out) :: distmin(nclosest)   ! the distance (m) of the closest column(s)
        1036             :     real(r8),intent(out) :: mlats(nclosest)     ! the latitude of the closest column(s)
        1037             :     real(r8),intent(out) :: mlons(nclosest)     ! the longitude of the closest column(s)
        1038             :
        1039             :     real(r8) :: dist
        1040             :     real(r8) :: rlats(pcols), rlons(pcols)
        1041             :     real(r8) :: latr, lonr
        1042             :
        1043           0 :     integer :: my_owner(nclosest)
        1044           0 :     integer :: my_lcid(nclosest)
        1045           0 :     integer :: my_icol(nclosest)
        1046           0 :     integer :: my_gcol(nclosest)
        1047           0 :     real(r8) :: my_distmin(nclosest)
        1048           0 :     real(r8) :: my_mlats(nclosest)
        1049           0 :     real(r8) :: my_mlons(nclosest)
        1050             :
        1051             :     integer  :: c, i, j, k, ierr, ncols, mindx(1)
        1052             :     real(r8) :: sendbufr(3)
        1053           0 :     real(r8) :: recvbufr(3,npes)
        1054             :     integer  :: sendbufi(4)
        1055           0 :     integer  :: recvbufi(4,npes)
        1056             :
        1057           0 :     call t_startf ('sat_hist::find_cols')
        1058             :
        1059           0 :     latr = lat/rad2deg              ! to radians
        1060           0 :     lonr = lon/rad2deg              ! to radians
        1061             :
        1062           0 :     my_owner(:)   = -999
        1063           0 :     my_lcid(:)    = -999
        1064           0 :     my_icol(:)    = -999
        1065           0 :     my_gcol(:)    = -999
        1066           0 :     my_mlats(:)   = -999
        1067           0 :     my_mlons(:)   = -999
        1068           0 :     my_distmin(:) = 1.e10_r8
        1069             :
        1070           0 :     chk_loop: do c=begchunk,endchunk
        1071           0 :        ncols = get_ncols_p(c)
        1072           0 :        call get_rlat_all_p(c, pcols, rlats)
        1073           0 :        call get_rlon_all_p(c, pcols, rlons)
        1074             :
        1075           0 :        col_loop: do i = 1,ncols
        1076             :           ! Use the Spherical Law of Cosines to find the great-circle distance.
        1077           0 :           dist = acos(sin(latr) * sin(rlats(i)) + cos(latr) * cos(rlats(i)) * cos(rlons(i) - lonr)) * rearth
        1078             :
        1079           0 :           closest_loop: do j = nclosest, 1, -1
        1080           0 :              if (dist < my_distmin(j)) then
        1081             :
        1082           0 :                 if (j < nclosest) then
        1083           0 :                    my_distmin(j+1) = my_distmin(j)
        1084           0 :                    my_owner(j+1)   = my_owner(j)
        1085           0 :                    my_lcid(j+1)    = my_lcid(j)
        1086           0 :                    my_icol(j+1)    = my_icol(j)
        1087           0 :                    my_gcol(j+1)    = my_gcol(j)
        1088           0 :                    my_mlats(j+1)   = my_mlats(j)
        1089           0 :                    my_mlons(j+1)   = my_mlons(j)
        1090             :                 end if
        1091             :
        1092           0 :                 my_distmin(j) = dist
        1093           0 :                 my_owner(j)   = iam
        1094           0 :                 my_lcid(j)    = c
        1095           0 :                 my_icol(j)    = i
        1096           0 :                 my_gcol(j)    = get_gcol_p(c,i)
        1097           0 :                 my_mlats(j)   = rlats(i) * rad2deg
        1098           0 :                 my_mlons(j)   = rlons(i) * rad2deg
        1099             :              else
        1100             :                 exit
        1101             :              end if
        1102             :           enddo closest_loop
        1103             :
        1104             :        enddo col_loop
        1105             :     enddo chk_loop
        1106             :
        1107             :     k = 1
        1108             :
        1109           0 :     do j = 1, nclosest
        1110             :
        1111           0 :        sendbufr(1) = my_distmin(k)
        1112           0 :        sendbufr(2) = my_mlats(k)
        1113           0 :        sendbufr(3) = my_mlons(k)
        1114             :
        1115           0 :        call mpi_allgather( sendbufr, 3, mpi_real8, recvbufr, 3, mpi_real8, mpicom, ierr )
        1116             :
        1117           0 :        mindx = minloc(recvbufr(1,:))
        1118           0 :        distmin(j) = recvbufr(1,mindx(1))
        1119           0 :        mlats(j)   = recvbufr(2,mindx(1))
        1120           0 :        mlons(j)   = recvbufr(3,mindx(1))
        1121             :
        1122           0 :        sendbufi(1) = my_owner(k)
        1123           0 :        sendbufi(2) = my_lcid(k)
        1124           0 :        sendbufi(3) = my_icol(k)
        1125           0 :        sendbufi(4) = my_gcol(k)
        1126             :
        1127           0 :        call mpi_allgather( sendbufi, 4, mpi_integer, recvbufi, 4, mpi_integer, mpicom, ierr )
        1128             :
        1129           0 :        owner(j)   = recvbufi(1,mindx(1))
        1130           0 :        lcid(j)    = recvbufi(2,mindx(1))
        1131           0 :        icol(j)    = recvbufi(3,mindx(1))
        1132           0 :        gcol(j)    = recvbufi(4,mindx(1))
        1133             :
        1134           0 :        if ( iam == owner(j) ) then
        1135           0 :           k = k+1
        1136             :        endif
        1137             :
        1138             :     enddo
        1139             :
        1140           0 :     call t_stopf ('sat_hist::find_cols')
        1141             :
        1142           0 :   end subroutine find_cols
        1143             :
        1144             : end module sat_hist

Generated by: LCOV version 1.14