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