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