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 32811448 : logical function is_satfile (file_index)
99 : integer, intent(in) :: file_index ! index of file in question
100 32811448 : is_satfile = file_index == sat_tape_num
101 32811448 : 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
|
