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