FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
interpolator.inc
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19!> @defgroup interpolator_mod interpolator_mod
20!> @ingroup interpolator
21!> @brief A module to interpolate climatology data to model the grid.
22!> @author William Cooke <William.Cooke@noaa.gov>
23!#######################################################################
24!
25!---------------------------------------------------------------------
26!> @brief interpolator_init receives various data as input in order
27!! to initialize interpolating.
28!!
29!! @param [inout] <clim_type> An interpolate type containing the necessary file and field
30!! data to be passed to the interpolator routine
31!! @param [in] <file_name> Climatology filename
32!! @param [in] <lonb_mod> The corners of the model grid-box longitudes
33!! @param [in] <latb_mod> The corners of the model grid_box latitudes
34!! @param [in] <data_names> OPTIONAL: A list of the names of components within the climatology
35!! file which you wish to read
36!! @param [in] <data_out_of_bounds> A list of the flags that are to be used in determining
37!! what to do if the pressure levels in the model go out of
38!! bounds from those of the climatology
39!! @param [in] <vert_interp> OPTIONAL: Flag to determine type of vertical interpolation
40!! @param [out] <clim_units> OPTIONAL: A list of the units for the components listed in data_names
41!! @param [out] <single_year_file> OPTIONAL: No description
42subroutine interpolator_init_( clim_type, file_name, lonb_mod, latb_mod, &
43 data_names, data_out_of_bounds, &
44 vert_interp, clim_units, single_year_file)
45type(interpolate_type), intent(inout) :: clim_type
46character(len=*), intent(in) :: file_name
47real(FMS_INTP_KIND_), intent(in) :: lonb_mod(:,:), latb_mod(:,:)
48character(len=*), intent(in) , optional :: data_names(:)
49!++lwh
50integer , intent(in) :: data_out_of_bounds(:)
51integer , intent(in), optional :: vert_interp(:)
52!--lwh
53character(len=*), intent(out), optional :: clim_units(:)
54logical, intent(out), optional :: single_year_file
55!
56! INTENT IN
57! file_name :: Climatology filename
58! lonb_mod :: The corners of the model grid-box longitudes.
59! latb_mod :: The corners of the model grid_box latitudes.
60! data_names :: A list of the names of components within the climatology file which you wish to read.
61! data_out_of_bounds :: A list of the flags that are to be used in determining what to do
62! if the pressure levels in the model
63! go out of bounds from those of the climatology.
64! vert_interp:: Flag to determine type of vertical interpolation
65!
66! INTENT OUT
67! clim_type :: An interpolate type containing the necessary file and field data to be passed
68! to the interpolator routine.
69! clim_units :: A list of the units for the components listed in data_names.
70!
71integer :: io, ierr
72
73if (.not. module_is_initialized) then
74 call fms_init
75 call diag_manager_init
76 call horiz_interp_init
77
78!--------------------------------------------------------------------
79! namelist input
80!--------------------------------------------------------------------
81
82 read (input_nml_file, nml=interpolator_nml, iostat=io)
83 ierr = check_nml_error(io,'interpolator_nml')
84
85 ! retain_cm3_bug is no longer supported.
86 if (retain_cm3_bug) then
87 call mpp_error(fatal, "interpolator_init: You have overridden the default value of " // &
88 "retain_cm3_bug and set it to .true. in interpolator_nml. This was a temporary workaround " // &
89 "that is no longer supported. Please remove this namelist variable.")
90 endif
91
92!---------------------------------------------------------------------
93! write version number and namelist to logfile.
94!---------------------------------------------------------------------
95call write_version_number("INTERPOLATOR_MOD", version)
96
97if (mpp_pe() == mpp_root_pe() ) write (stdlog(), nml=interpolator_nml)
98
99module_is_initialized = .true.
100
101endif !> if (module_is_initilized)
102if (use_mpp_io) then
103 call mpp_error(fatal, "Interpolator::nml=interpolator_nml " //&
104 'MPP_IO is no longer supported. Please remove from use_mpp_io from interpolator_nml')
105endif
106
107clim_type%FMS_INTP_TYPE_%is_allocated = .true.
108
109call fms2io_interpolator_init(clim_type, file_name, lonb_mod, latb_mod, &
110 data_names, data_out_of_bounds, &
111 vert_interp, clim_units, single_year_file)
112
113end subroutine interpolator_init_
114
115subroutine fms2io_interpolator_init_(clim_type, file_name, lonb_mod, latb_mod, &
116 data_names, data_out_of_bounds, &
117 vert_interp, clim_units, single_year_file)
118
119type(interpolate_type), intent(inout) :: clim_type
120character(len=*), intent(in) :: file_name
121real(FMS_INTP_KIND_) , intent(in) :: lonb_mod(:,:), latb_mod(:,:)
122character(len=*), intent(in) , optional :: data_names(:)
123!++lwh
124integer , intent(in) :: data_out_of_bounds(:)
125integer , intent(in), optional :: vert_interp(:)
126!--lwh
127character(len=*), intent(out), optional :: clim_units(:)
128logical, intent(out), optional :: single_year_file
129
130character(len=FMS_FILE_LEN) :: src_file
131!++lwh
132real(FMS_INTP_KIND_) :: dlat, dlon
133!--lwh
134type(time_type) :: base_time
135integer :: fileday, filemon, fileyr, filehr, filemin,filesec, m,m1
136character(len= 20) :: fileunits
137character(len=128) :: var_dimname(6)
138character(len=128), allocatable :: var_names(:)
139integer, allocatable :: var_ndims(:)
140integer :: j, i
141logical :: non_monthly
142character(len=24) :: file_calendar
143character(len=256) :: error_mesg
144integer :: model_calendar
145integer :: yr, mo, dy, hr, mn, sc
146integer :: n
147type(time_type) :: Julian_time, Noleap_time
148real(r8_kind), allocatable :: time_in(:)
149real(FMS_INTP_KIND_), allocatable, save :: agrid_mod(:,:,:)
150integer :: nx, ny
151type(FmsNetcdfFile_t) :: fileobj
152
153integer, parameter :: lkind=fms_intp_kind_
154real(FMS_INTP_KIND_), parameter :: lPI=real(pi,fms_intp_kind_)
155
156!> variables used to set time
157logical :: yearly !< flags to indicate if time data is in units of months or years
158integer :: num_years !< number of years
159integer :: base_year, base_month, base_day, base_hour, base_minute, base_second
160integer :: nn !< counter
161logical :: noleap_file_calendar !< is the file calendar noleap or julian
162real(r8_kind) :: num_days, frac_year !< variables for yearly=.true.
163type(time_type) :: n_time !< temporary time
164
165clim_type%separate_time_vary_calc = .false.
166
167num_fields = 0
168
169!--------------------------------------------------------------------
170! open source file containing fields to be interpolated
171!--------------------------------------------------------------------
172src_file = 'INPUT/'//trim(file_name)
173
174if(fms2_io_file_exist(trim(src_file))) then
175 if(.not. open_file(clim_type%fileobj, trim(src_file), 'read')) &
176 call mpp_error(fatal, 'Interpolator_init: Error in opening file '//trim(src_file))
177 fileobj = clim_type%fileobj
178else
179!Climatology file doesn't exist, so exit
180 call mpp_error(fatal,'Interpolator_init : Data file '//trim(src_file)//' does not exist')
181endif
182
183!Find the number of variables (nvar) in this file
184clim_type%file_name = trim(file_name)
185nvar = get_num_variables(fileobj)
186num_fields = nvar
187if(present(data_names)) then
188 num_fields= size(data_names(:))
189else
190 allocate(var_names(nvar), var_ndims(nvar))
191 call get_variable_names(fileobj, var_names)
192 !--- loop through all the vars to exclude scalar or 1-D array
193 do i=1,nvar
194 var_ndims(i) = get_variable_num_dimensions(fileobj, var_names(i))
195 enddo
196 num_fields = count(var_ndims>1)
197endif
198
199! -------------------------------------------------------------------
200! Allocate space for the number of axes in the data file.
201! -------------------------------------------------------------------
202
203nlon=0 ! Number of longitudes (center-points) in the climatology.
204nlat=0 ! Number of latitudes (center-points) in the climatology.
205nlev=0 ! Number of levels (center-points) in the climatology.
206nlatb=0 ! Number of longitudes (boundaries) in the climatology.
207nlonb=0 ! Number of latitudes (boundaries) in the climatology.
208nlevh=0 ! Number of levels (boundaries) in the climatology.
209
210clim_type%level_type = 0 ! Default value
211
212!++lwh
213! -------------------------------------------------------------------
214! For 2-D fields, set a default value of nlev=nlevh=1
215! -------------------------------------------------------------------
216nlev = 1
217nlevh = 1
218!--lwh
219 clim_type%vertical_indices = 0 ! initial value
220
221!--- get clim_type%lat
222if(dimension_exists(fileobj, "lat")) then
223 call get_dimension_size(fileobj, "lat", nlat)
224else
225 call mpp_error(fatal,'Interpolator_init : dimension lat does not exist in file '//trim(src_file) )
226endif
227allocate(clim_type%FMS_INTP_TYPE_%lat(nlat))
228call get_axis_latlon_data(fileobj, 'lat', clim_type%FMS_INTP_TYPE_%lat)
229
230!--- get clim_type%lon
231if(dimension_exists(fileobj, "lon")) then
232 call get_dimension_size(fileobj, "lon", nlon)
233else
234 call mpp_error(fatal,'Interpolator_init : dimension lon does not exist in file '//trim(src_file) )
235endif
236allocate(clim_type%FMS_INTP_TYPE_%lon(nlon))
237call get_axis_latlon_data(fileobj, 'lon', clim_type%FMS_INTP_TYPE_%lon)
238
239!--- get clim_type%latb
240if(dimension_exists(fileobj, "latb")) then
241 call get_dimension_size(fileobj, "latb", nlatb)
242 allocate(clim_type%FMS_INTP_TYPE_%latb(nlatb))
243 call get_axis_latlon_data(fileobj, 'latb', clim_type%FMS_INTP_TYPE_%latb)
244else
245 if(nlat == 1) call mpp_error(fatal,'Interpolator_init : nlat is 1')
246 ! In the case where only the grid midpoints of the latitudes are defined we force the
247 ! definition of the boundaries to be half-way between the midpoints.
248 allocate(clim_type%FMS_INTP_TYPE_%latb(nlat+1))
249 dlat = (clim_type%FMS_INTP_TYPE_%lat(2)-clim_type%FMS_INTP_TYPE_%lat(1)) * 0.5_lkind
250! clim_type%latb(1) = min( 90., max(-90., clim_type%lat(1) - dlat) )
251 clim_type%FMS_INTP_TYPE_%latb(1) = min( lpi/2._lkind, max(-lpi/2._lkind, clim_type%FMS_INTP_TYPE_%lat(1) - dlat) )
252 clim_type%FMS_INTP_TYPE_%latb(2:nlat) = ( clim_type%FMS_INTP_TYPE_%lat(1:nlat-1) + &
253 clim_type%FMS_INTP_TYPE_%lat(2:nlat) ) * 0.5_lkind
254 dlat = ( clim_type%FMS_INTP_TYPE_%lat(nlat) - clim_type%FMS_INTP_TYPE_%lat(nlat-1) ) * 0.5_lkind
255! clim_type%latb(nlat+1) = min( 90., max(-90., clim_type%lat(nlat) + dlat) )
256 clim_type%FMS_INTP_TYPE_%latb(nlat+1) = min( lpi/2._lkind, max(-lpi/2._lkind, &
257 clim_type%FMS_INTP_TYPE_%lat(nlat) + dlat) )
258endif
259
260!--- get clim_type%lonb
261if(dimension_exists(fileobj, "lonb")) then
262 call get_dimension_size(fileobj, "lonb", nlonb)
263 allocate(clim_type%FMS_INTP_TYPE_%lonb(nlonb))
264 call get_axis_latlon_data(fileobj, 'lonb', clim_type%FMS_INTP_TYPE_%lonb)
265else
266! In the case where only the midpoints of the longitudes are defined we force the definition
267! of the boundaries to be half-way between the midpoints.
268 if (size(clim_type%FMS_INTP_TYPE_%lon(:)) /= 1) then
269 allocate(clim_type%FMS_INTP_TYPE_%lonb(size(clim_type%FMS_INTP_TYPE_%lon(:))+1))
270 dlon = ( clim_type%FMS_INTP_TYPE_%lon(2)-clim_type%FMS_INTP_TYPE_%lon(1) )/2.0_lkind
271 clim_type%FMS_INTP_TYPE_%lonb(1) = clim_type%FMS_INTP_TYPE_%lon(1) - dlon
272 clim_type%FMS_INTP_TYPE_%lonb(2:) = clim_type%FMS_INTP_TYPE_%lon(1:) + dlon
273 else
274!! this is the case for zonal mean data, lon = 1, lonb not present
275!! in file.
276 allocate(clim_type%FMS_INTP_TYPE_%lonb(2))
277 clim_type%FMS_INTP_TYPE_%lonb(1) = -360.0_lkind*real(dtr,fms_intp_kind_)
278 clim_type%FMS_INTP_TYPE_%lonb(2) = 360.0_lkind*real(dtr,fms_intp_kind_)
279 clim_type%FMS_INTP_TYPE_%lon(1) = 0.0_lkind
280 endif
281endif
282
283!--- get clim_type%levs
284clim_type%level_type = 0
285clim_type%vertical_indices = 0
286if(dimension_exists(fileobj, "pfull")) then
287 call get_dimension_size(fileobj, "pfull", nlev)
288 allocate(clim_type%FMS_INTP_TYPE_%levs(nlev))
289 call get_axis_level_data(fileobj, 'pfull', clim_type%FMS_INTP_TYPE_%levs, &
290 clim_type%level_type, clim_type%vertical_indices)
291else if(dimension_exists(fileobj, "sigma_full")) then
292 call get_dimension_size(fileobj, "sigma_full", nlev)
293 allocate(clim_type%FMS_INTP_TYPE_%levs(nlev))
294 call fms2_io_read_data(fileobj, "sigma_full", clim_type%FMS_INTP_TYPE_%levs)
295 clim_type%level_type = sigma
296else
297 ! -------------------------------------------------------------------
298 ! For 2-D fields, allocate levs and halflevs here
299 ! code is still needed for case when only halflevs are in data file.
300 ! -------------------------------------------------------------------
301 nlev = 1
302 allocate( clim_type%FMS_INTP_TYPE_%levs(nlev) )
303 clim_type%FMS_INTP_TYPE_%levs = 0.0_lkind
304endif
305
306!--- get clim_type%halflevs
307if(dimension_exists(fileobj, "phalf")) then
308 call get_dimension_size(fileobj, "phalf", nlevh)
309 allocate(clim_type%FMS_INTP_TYPE_%halflevs(nlevh))
310 call get_axis_level_data(fileobj, 'phalf', clim_type%FMS_INTP_TYPE_%halflevs, &
311 clim_type%level_type, clim_type%vertical_indices)
312else if(dimension_exists(fileobj, "sigma_half")) then
313 call get_dimension_size(fileobj, "sigma_half", nlevh)
314 allocate(clim_type%FMS_INTP_TYPE_%halflevs(nlevh))
315 call fms2_io_read_data(fileobj, "sigma_half", clim_type%FMS_INTP_TYPE_%halflevs)
316 clim_type%level_type = sigma
317else
318 allocate( clim_type%FMS_INTP_TYPE_%halflevs(nlev+1) )
319 clim_type%FMS_INTP_TYPE_%halflevs(1) = 0.0_lkind
320 if (clim_type%level_type == pressure) then
321 clim_type%FMS_INTP_TYPE_%halflevs(nlev+1) = 1013.25_lkind* 100.0_lkind ! MKS
322 else if (clim_type%level_type == sigma ) then
323 clim_type%FMS_INTP_TYPE_%halflevs(nlev+1) = 1.0_lkind
324 endif
325 do n=2,nlev
326 clim_type%FMS_INTP_TYPE_%halflevs(n) = 0.5_lkind*(clim_type%FMS_INTP_TYPE_%levs(n) + &
327 clim_type%FMS_INTP_TYPE_%levs(n-1))
328 end do
329endif
330
331!get time informaiton
332if(dimension_exists(fileobj, "time")) then
333 call get_dimension_size(fileobj, "time", ntime)
334
335 call get_variable_units(fileobj, "time", units)
336 call get_time_calendar(fileobj, "time", file_calendar)
337 model_calendar = get_calendar_type()
338 fileday = 0
339 filemon = 0
340 fileyr = 0
341 filehr = 0
342 filemin = 0
343 filesec = 0
344 yearly = .false.
345 select case(units(:3))
346 case('day')
347 fileunits = units(12:) !Assuming "days since YYYY-MM-DD HH:MM:SS"
348 if ( len_trim(fileunits) < 19 ) then
349 write(error_mesg, '(A49,A,A49,A)' ) &
350 'Interpolator_init : Incorrect time units in file ', &
351 trim(file_name), '. Expecting days since YYYY-MM-DD HH:MM:SS, found', &
352 trim(units)
353 call mpp_error(fatal,error_mesg)
354 endif
355 read(fileunits(1:4) , *) fileyr
356 read(fileunits(6:7) , *) filemon
357 read(fileunits(9:10) , *) fileday
358 read(fileunits(12:13), *) filehr
359 read(fileunits(15:16), *) filemin
360 read(fileunits(18:19), *) filesec
361 case('mon')
362 fileunits = units(14:) !Assuming "months since YYYY-MM-DD HH:MM:SS"
363 if ( len_trim(fileunits) < 19 ) then
364 write(error_mesg, '(A49,A,A51,A)' ) &
365 'Interpolator_init : Incorrect time units in file ', &
366 trim(file_name), '. Expecting months since YYYY-MM-DD HH:MM:SS, found', &
367 trim(units)
368 call mpp_error(fatal,error_mesg)
369 endif
370 read(fileunits(1:4) , *) fileyr
371 read(fileunits(6:7) , *) filemon
372 read(fileunits(9:10) , *) fileday
373 read(fileunits(12:13), *) filehr
374 read(fileunits(15:16), *) filemin
375 read(fileunits(18:19), *) filesec
376 case( 'yea')
377 fileunits = units(13:) !Assuming "years since YYYY-MM-DD HH:MM:SS"
378 if ( len_trim(fileunits) < 19 ) then
379 write(error_mesg, '(A49,A,A51,A)' ) &
380 'Interpolator_init : Incorrect time units in file ', &
381 trim(file_name), '. Expecting years since YYYY-MM-DD HH:MM:SS, found', &
382 trim(units)
383 call mpp_error(fatal,error_mesg)
384 endif
385 read(fileunits(1:4) , *) fileyr
386 read(fileunits(6:7) , *) filemon
387 read(fileunits(9:10) , *) fileday
388 read(fileunits(12:13), *) filehr
389 read(fileunits(15:16), *) filemin
390 read(fileunits(18:19), *) filesec
391 yearly = .true.
392 case default
393 call mpp_error(fatal,'Interpolator_init : Time units not recognized in file '//file_name)
394 end select
395
396 clim_type%climatological_year = (fileyr == 0)
397
398 if (.not. clim_type%climatological_year) then
399
400 !----------------------------------------------------------------------
401 ! if file date has a non-zero year in the base time, determine that
402 ! base_time based on the netcdf info.
403 !----------------------------------------------------------------------
404 noleap_file_calendar=.false.
405 if ( (model_calendar == julian .and. trim(adjustl(lowercase(file_calendar))) == 'julian')) then
406 call mpp_error (note, 'interpolator_mod: Model and file&
407 & calendars are the same for file ' // &
408 & trim(file_name) // '; no calendar conversion &
409 &needed')
410 base_time = set_date(fileyr, filemon, fileday, filehr,filemin,filesec)
411 noleap_file_calendar=.false.
412 else if((model_calendar == noleap .and. trim(adjustl(lowercase(file_calendar))) == 'noleap') ) then
413 call mpp_error (note, 'interpolator_mod: Model and file&
414 & calendars are the same for file ' // &
415 & trim(file_name) // '; no calendar conversion &
416 &needed')
417 base_time = set_date(fileyr, filemon, fileday, filehr,filemin,filesec)
418 noleap_file_calendar=.true.
419 else if ( (model_calendar == julian .and. trim(adjustl(lowercase(file_calendar))) == 'noleap')) then
420 call mpp_error (note, 'interpolator_mod: Using julian &
421 &model calendar and noleap file calendar&
422 & for file ' // trim(file_name) // &
423 &'; calendar conversion needed')
424 base_time = set_date_no_leap(fileyr, filemon, fileday,filehr, filemin, filesec)
425 noleap_file_calendar=.true.
426 else if ( (model_calendar == noleap .and. trim(adjustl(lowercase(file_calendar))) == 'julian')) then
427 call mpp_error (note, 'interpolator_mod: Using noleap &
428 &model calendar and julian file calendar&
429 & for file ' // trim(file_name) // &
430 &'; calendar conversion needed')
431 base_time = set_date_julian(fileyr, filemon, fileday,filehr, filemin, filesec)
432 noleap_file_calendar=.false.
433 else
434 call mpp_error (fatal , 'interpolator_mod: Model and file&
435 & calendars ( ' // trim(file_calendar) // ' ) differ &
436 &for file ' // trim(file_name) // '; this calendar &
437 &conversion not currently available')
438 endif
439
440 else
441
442 !! if the year is specified as '0000', then the file is intended to
443 !! apply to all years -- the time variables within the file refer to
444 !! the displacement from the start of each year to the time of the
445 !! associated data. Time interpolation is to be done with interface
446 !! time_interp_list, with the optional argument modtime=YEAR. base_time
447 !! is set to an arbitrary value here; it's only use will be as a
448 !! timestamp for optionally generated diagnostics.
449 base_time = get_base_time()
450 endif
451
452
453 ntime_in = 1
454 if (ntime > 0) then
455 allocate(time_in(ntime), clim_type%time_slice(ntime))
456 allocate(clim_type%clim_times(12,(ntime+11)/12))
457 time_in = 0.0_r8_kind
458 clim_type%time_slice = set_time(0,0) + base_time
459 clim_type%clim_times = set_time(0,0) + base_time
460 call fms2_io_read_data(fileobj, "time", time_in)
461 ntime_in = ntime
462 !> convert the number of years passed to days if yearly=.true.
463 if(yearly) then
464 do n = 1, ntime
465 if(.not.noleap_file_calendar) then
466 ! Julian calendar
467 num_years = int(time_in(n))
468 frac_year = time_in(n) - real(num_years, r8_kind)
469 call get_date_julian(base_time, base_year, base_month, base_day, base_hour, base_minute, base_second)
470 num_days = 0.0_r8_kind
471 do nn=1, num_years
472 if( mod(base_year+nn-1,4)==0) num_days = num_days + 366._r8_kind
473 if( mod(base_year+nn-1,4)/=0) num_days = num_days + 365._r8_kind
474 end do
475 if( mod(base_year+num_years,4)==0) num_days = num_days + 366._r8_kind*frac_year
476 if( mod(base_year+num_years,4)/=0) num_days = num_days + 365._r8_kind*frac_year
477 else
478 num_days = time_in(n)*365._r8_kind
479 end if
480 time_in(n)=num_days
481 end do
482 end if
483 ! determine whether the data is a continuous set of monthly values or
484 ! a series of annual cycles spread throughout the period of data
485 non_monthly = .false.
486 do n = 1, ntime-1
487 ! Assume that the times in the data file correspond to days only.
488 if (time_in(n+1) > (time_in(n) + 32._r8_kind)) then
489 non_monthly = .true.
490 exit
491 endif
492 end do
493 if (clim_type%climatological_year) then
494 call mpp_error (note, 'interpolator_mod :' // &
495 & trim(file_name) // ' is a year-independent climatology file')
496 else
497 call mpp_error (note, 'interpolator_mod :' // &
498 & trim(file_name) // ' is a timeseries file')
499 endif
500
501 do n = 1, ntime
502 !Assume that the times in the data file correspond to days only.
503
504
505 if (clim_type%climatological_year) then
506 !! RSH NOTE:
507 !! for this case, do not add base_time. time_slice will be sent to
508 !! time_interp_list with the optional argument modtime=YEAR, so that
509 !! the time that is needed in time_slice is the displacement into the
510 !! year, not the displacement from a base_time.
511 clim_type%time_slice(n) = set_time( int( (time_in(n)-real(int(time_in(n)),r8_kind))*seconds_per_day), &
512 int(time_in(n)))
513 else
514
515 !--------------------------------------------------------------------
516 ! if fileyr /= 0 (i.e., climatological_year=F),
517 ! then define the times associated with each time-
518 ! slice. if calendar conversion between data file and model calendar
519 ! is needed, do it so that data from the file is associated with the
520 ! same calendar time in the model. here the time_slice needs to
521 ! include the base_time; values will be generated relative to the
522 ! "real" time.
523 !--------------------------------------------------------------------
524 if ( (model_calendar == julian .and. trim(adjustl(lowercase(file_calendar))) == 'julian') .or. &
525 & (model_calendar == noleap .and. trim(adjustl(lowercase(file_calendar))) == 'noleap') ) then
526
527 !---------------------------------------------------------------------
528 ! no calendar conversion needed.
529 !---------------------------------------------------------------------
530 clim_type%time_slice(n) = &
531 set_time( int( (time_in(n)-real(int(time_in(n)),r8_kind))*seconds_per_day), &
532 int(time_in(n))) + base_time
533
534 !---------------------------------------------------------------------
535 ! convert file times from noleap to julian.
536 !---------------------------------------------------------------------
537 else if ( (model_calendar == julian .and. trim(adjustl(lowercase(file_calendar))) == 'noleap')) then
538 noleap_time = set_time(0, int(time_in(n))) + base_time
539 call get_date_no_leap (noleap_time, yr, mo, dy, hr, mn, sc)
540 clim_type%time_slice(n) = set_date_julian(yr, mo, dy, hr, mn, sc)
541 if (n == 1) then
542 call print_date (clim_type%time_slice(1), &
543 str= 'for file ' // trim(file_name) // ', the first time slice is mapped to :')
544 endif
545 if (n == ntime) then
546 call print_date (clim_type%time_slice(ntime), &
547 str= 'for file ' // trim(file_name) // ', the last time slice is mapped to:')
548 endif
549
550
551 !---------------------------------------------------------------------
552 ! convert file times from julian to noleap.
553 !---------------------------------------------------------------------
554 else if ( (model_calendar == noleap .and. trim(adjustl(lowercase(file_calendar))) == 'julian')) then
555 julian_time = set_time(0, int(time_in(n))) + base_time
556 call get_date_julian (julian_time, yr, mo, dy, hr, mn, sc)
557 clim_type%time_slice(n) = set_date_no_leap(yr, mo, dy,hr, mn, sc)
558 if (n == 1) then
559 call print_date (clim_type%time_slice(1), &
560 str= 'for file ' // trim(file_name) // ', the first time slice is mapped to :')
561 endif
562 if (n == ntime) then
563 call print_date (clim_type%time_slice(ntime), &
564 str= 'for file ' // trim(file_name) // ', the last time slice is mapped to:')
565 endif
566
567 !---------------------------------------------------------------------
568 ! any other calendar combinations would have caused a fatal error
569 ! above.
570 !---------------------------------------------------------------------
571 endif
572 endif
573
574 m = (n-1)/12 +1 ; m1 = n- (m-1)*12
575 clim_type%clim_times(m1,m) = clim_type%time_slice(n)
576 enddo
577 else
578 allocate(time_in(1), clim_type%time_slice(1))
579 allocate(clim_type%clim_times(1,1))
580 time_in = 0.0_r8_kind
581 clim_type%time_slice = set_time(0,0) + base_time
582 clim_type%clim_times(1,1) = set_time(0,0) + base_time
583 endif
584 deallocate(time_in)
585else
586 call mpp_error(fatal, 'Interpolator_init: time axis does not exist in file '//trim(src_file))
587endif
588
589!Assume that the horizontal interpolation within a file is the same for each variable.
590
591 if (conservative_interp) then
592 call horiz_interp_new (clim_type%interph, &
593 clim_type%FMS_INTP_TYPE_%lonb, clim_type%FMS_INTP_TYPE_%latb, &
594 lonb_mod, latb_mod)
595 else
596
597 call mpp_error(note, "Using Bilinear interpolation")
598
599 !!! DEBUG CODE
600 if (.not. allocated(agrid_mod)) then
601 nx = size(lonb_mod,1)-1
602 ny = size(latb_mod,2)-1
603 allocate(agrid_mod(nx,ny,2))
604 do j=1,ny
605 do i=1,nx
606 call cell_center2((/lonb_mod(i,j),latb_mod(i,j)/), &
607 (/lonb_mod(i+1,j),latb_mod(i+1,j)/), &
608 (/lonb_mod(i,j+1),latb_mod(i,j+1)/), &
609 (/lonb_mod(i+1,j+1),latb_mod(i+1,j+1)/), agrid_mod(i,j,:))
610 enddo
611 enddo
612 endif
613
614 !!! END DEBUG CODE
615
616 call horiz_interp_new (clim_type%interph, &
617 clim_type%FMS_INTP_TYPE_%lonb, clim_type%FMS_INTP_TYPE_%latb, &
618 agrid_mod(:,:,1), agrid_mod(:,:,2), interp_method="bilinear")
619 endif
620
621!--------------------------------------------------------------------
622! allocate the variable clim_type%data . This will be the climatology
623! data horizontally interpolated, so it will be on the model horizontal
624! grid, but it will still be on the climatology vertical grid.
625!--------------------------------------------------------------------
626
627select case(ntime)
628 case (13:)
629! This may be data that does not have a continous time-line
630! i.e. IPCC data where decadal data is present but we wish to retain
631! the seasonal nature of the data.
632!! RSH: the following test will not always work; instead use the
633!! RSH: non-monthly variable to test on.
634!RSHlast_time = clim_type%time_slice(1) + ( ntime -1 ) * &
635!RSH ( clim_type%time_slice(2) - clim_type%time_slice(1) )
636
637!RSHif ( last_time < clim_type%time_slice(ntime)) then
638
639 if (non_monthly) then
640! We have a broken time-line. e.g. We have monthly data but only for years ending in 0. 1960,1970 etc.
641! allocate(clim_type%data5d(size(lonb_mod(:))-1, size(latb_mod(:))-1, nlev, 2, num_fields))
642 allocate(clim_type%FMS_INTP_TYPE_%pmon_pyear(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, num_fields))
643 allocate(clim_type%FMS_INTP_TYPE_%pmon_nyear(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, num_fields))
644 allocate(clim_type%FMS_INTP_TYPE_%nmon_nyear(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, num_fields))
645 allocate(clim_type%FMS_INTP_TYPE_%nmon_pyear(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, num_fields))
646 clim_type%FMS_INTP_TYPE_%pmon_pyear = 0.0_lkind
647 clim_type%FMS_INTP_TYPE_%pmon_nyear = 0.0_lkind
648 clim_type%FMS_INTP_TYPE_%nmon_nyear = 0.0_lkind
649 clim_type%FMS_INTP_TYPE_%nmon_pyear = 0.0_lkind
650 clim_type%TIME_FLAG = bilinear
651else
652! We have a continuous time-line so treat as for 5-12 timelevels as below.
653 if ( .not. read_all_on_init) then
654 allocate(clim_type%FMS_INTP_TYPE_%data5d(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, 2, num_fields))
655 else
656 allocate(clim_type%FMS_INTP_TYPE_%data5d(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, &
657 ntime, num_fields))
658 endif
659 clim_type%FMS_INTP_TYPE_%data5d = 0.0_lkind
660 clim_type%TIME_FLAG = linear
661endif
662
663
664!++lwh
665 case (1:12)
666!--lwh
667! We have more than 4 timelevels
668! Assume we have monthly or higher time resolution datasets (climatology or time series)
669! So we only need to read 2 datasets and apply linear temporal interpolation.
670 if ( .not. read_all_on_init) then
671 allocate(clim_type%FMS_INTP_TYPE_%data5d(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, 2, num_fields))
672 else
673 allocate(clim_type%FMS_INTP_TYPE_%data5d(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, &
674 ntime, num_fields))
675 endif
676 clim_type%FMS_INTP_TYPE_%data5d = 0.0_lkind
677 clim_type%TIME_FLAG = linear
678!++lwh
679!case (1:4)
680! Assume we have seasonal data and read in all the data.
681! We can apply sine curves to these data.
682
683! allocate(clim_type%data5d(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, ntime, num_fields))
684! clim_type%data5d = 0.0
685! clim_type%TIME_FLAG = SEASONAL
686!--lwh
687! case (default)
688 case(:0)
689 clim_type%TIME_FLAG = notime
690 allocate(clim_type%FMS_INTP_TYPE_%data5d(size(lonb_mod,1)-1, size(latb_mod,2)-1, nlev, 1, num_fields))
691end select
692
693
694!------------------------------------------------------------------
695! Allocate space for the single time level of the climatology on its
696! grid size.
697!----------------------------------------------------------------------
698
699 if(clim_type%TIME_FLAG .eq. linear ) then
700 allocate(clim_type%time_init(num_fields,2))
701 else
702 allocate(clim_type%time_init(num_fields,ntime))
703 endif
704 allocate (clim_type%indexm(num_fields), &
705 clim_type%indexp(num_fields), &
706 clim_type%climatology(num_fields))
707 clim_type%time_init(:,:) = 0
708 clim_type%indexm(:) = 0
709 clim_type%indexp(:) = 0
710 clim_type%climatology(:) = 0
711
712allocate(clim_type%has_level(num_fields))
713allocate(clim_type%field_name(num_fields))
714allocate(clim_type%mr(num_fields))
715allocate(clim_type%out_of_bounds(num_fields))
716clim_type%out_of_bounds(:)=0
717allocate(clim_type%vert_interp(num_fields))
718clim_type%vert_interp(:)=0
719
720if(present(data_names)) then
721
722!++lwh
723 if ( size(data_out_of_bounds(:)) /= size(data_names(:)) .and. size(data_out_of_bounds(:)) /= 1 ) &
724 call mpp_error(fatal,'interpolator_init : The size of the data_out_of_bounds array must be 1&
725 & or size(data_names)')
726 if (present(vert_interp)) then
727 if( size(vert_interp(:)) /= size(data_names(:)) .and. size(vert_interp(:)) /= 1 ) &
728 call mpp_error(fatal,'interpolator_init : The size of the vert_interp array must be 1&
729 & or size(data_names)')
730 endif
731! Only read the fields named in data_names
732 do j=1,size(data_names(:))
733 if(variable_exists(fileobj, data_names(j)) ) then
734 call get_variable_units(fileobj, data_names(j), units)
735 ndim = get_variable_num_dimensions(fileobj, data_names(j))
736 clim_type%has_level(j) = .false.
737 if(ndim > 2) then
738 call get_variable_dimension_names(fileobj, data_names(j), var_dimname(1:ndim))
739 if(trim(var_dimname(3)) == "pfull" .OR. trim(var_dimname(3)) == "sigma_full") &
740 & clim_type%has_level(j) = .true.
741 endif
742
743 units=chomp(units)
744 if (mpp_pe() == 0 ) write(*,*) 'Initializing src field : ',trim(data_names(j))
745 clim_type%field_name(j) = data_names(j)
746 clim_type%mr(j) = check_climo_units(units)
747 if (present(clim_units)) clim_units(j) = units
748 clim_type%out_of_bounds(j) = data_out_of_bounds( min(j,SIZE(data_out_of_bounds(:))) )
749 if( clim_type%out_of_bounds(j) /= constant .and. clim_type%out_of_bounds(j) /= zero ) &
750 call mpp_error(fatal,"Interpolator_init: data_out_of_bounds must be&
751 & set to ZERO or CONSTANT")
752 if( present(vert_interp) ) then
753 clim_type%vert_interp(j) = vert_interp( min(j,SIZE(vert_interp(:))) )
754 if( clim_type%vert_interp(j) /= interp_weighted_p .and. clim_type%vert_interp(j) /= interp_linear_p ) &
755 call mpp_error(fatal,"Interpolator_init: vert_interp must be&
756 & set to INTERP_WEIGHTED_P or INTERP_LINEAR_P")
757 else
758 clim_type%vert_interp(j) = interp_weighted_p
759 end if
760 else
761 call mpp_error(fatal,'interpolator_init : Check names of fields being passed. ' &
762 &//trim(data_names(j))//' does not exist.')
763 endif
764 enddo
765else
766 if ( size(data_out_of_bounds(:)) /= nvar .and. size(data_out_of_bounds(:)) /= 1 ) &
767 call mpp_error(fatal,'interpolator_init : The size of the out of bounds array must be 1&
768 & or the number of fields in the climatology dataset')
769 if ( present(vert_interp) ) then
770 if (size(vert_interp(:)) /= nvar .and. size(vert_interp(:)) /= 1 ) &
771 call mpp_error(fatal,'interpolator_init : The size of the vert_interp array must be 1&
772 & or the number of fields in the climatology dataset')
773 endif
774
775! Read all the fields within the climatology data file.
776 j = 0
777 do i=1,nvar
778 ndim = var_ndims(i)
779 if(ndim .LE. 1) cycle
780 j = j + 1
781 clim_type%has_level(j) = .false.
782 if(ndim > 2) then
783 call get_variable_dimension_names(fileobj, var_names(i), var_dimname(1:ndim))
784 if(trim(var_dimname(3)) == "pfull" .OR. trim(var_dimname(3)) == "sigma_full") clim_type%has_level(j) = .true.
785 endif
786
787 call get_variable_units(fileobj, var_names(i), units)
788 if (mpp_pe() ==0 ) write(*,*) 'Initializing src field : ',trim(var_names(i))
789 clim_type%field_name(j) = trim(var_names(i))
790 clim_type%mr(j) = check_climo_units(units)
791 if (present(clim_units)) clim_units(j) = units
792 clim_type%out_of_bounds(j) = data_out_of_bounds( min(i,SIZE(data_out_of_bounds(:))) )
793 if( clim_type%out_of_bounds(j) /= constant .and. &
794 clim_type%out_of_bounds(j) /= zero ) &
795 call mpp_error(fatal,"Interpolator_init: data_out_of_bounds must be&
796 & set to ZERO or CONSTANT")
797 if( present(vert_interp) ) then
798 clim_type%vert_interp(j) = vert_interp( min(i,SIZE(vert_interp(:))) )
799 if( clim_type%vert_interp(j) /= interp_weighted_p .and. &
800 clim_type%vert_interp(j) /= interp_linear_p ) &
801 call mpp_error(fatal,"Interpolator_init: vert_interp must be&
802 & set to INTERP_WEIGHTED_P or INTERP_LINEAR_P")
803 else
804 clim_type%vert_interp(j) = interp_weighted_p
805 end if
806 end do
807 if(j .NE. num_fields) call mpp_error(fatal,"Interpolator_init: j does not equal to num_fields")
808 deallocate(var_names, var_ndims)
809!--lwh
810endif
811
812if( clim_type%TIME_FLAG .eq. seasonal ) then
813! Read all the data at this point.
814 do i=1,num_fields
815 do n = 1, ntime
816 call read_data( clim_type, clim_type%field_name(i), &
817 clim_type%FMS_INTP_TYPE_%data5d(:,:,:,n,i), n, i, base_time )
818 enddo
819 enddo
820endif
821
822if( clim_type%TIME_FLAG .eq. linear .and. read_all_on_init) then
823! Read all the data at this point.
824 do i=1,num_fields
825 do n = 1, ntime
826 call read_data( clim_type, clim_type%field_name(i), &
827 clim_type%FMS_INTP_TYPE_%data5d(:,:,:,n,i), n, i, base_time )
828 enddo
829 enddo
830
831 call close_file (fileobj)
832endif
833
834if( clim_type%TIME_FLAG .eq. notime ) then
835! Read all the data at this point.
836 do i=1,num_fields
837 call read_data_no_time_axis( clim_type, clim_type%field_name(i), &
838 clim_type%FMS_INTP_TYPE_%data5d(:,:,:,1,i), i )
839 enddo
840 call close_file (fileobj)
841endif
842
843if (present (single_year_file)) then
844 single_year_file = clim_type%climatological_year
845endif
846
847end subroutine fms2io_interpolator_init_
848
849subroutine get_axis_latlon_data_(fileobj, name, latlon_data)
850 type(FmsNetcdfFile_t), intent(in) :: fileobj
851 character(len=*), intent(in) :: name
852 real(FMS_INTP_KIND_), dimension(:), intent(out) :: latlon_data
853
854
855 if(variable_exists(fileobj, name)) then
856 call fms2_io_read_data(fileobj, name, latlon_data)
857 else
858 call mpp_error(fatal,'get_axis_latlon_data: variable '// &
859 & trim(name)//' does not exist in file '//trim(fileobj%path) )
860 endif
861 call get_variable_units(fileobj, name, units)
862 select case(units(1:6))
863 case('degree')
864 latlon_data = latlon_data*real(dtr,fms_intp_kind_)
865 case('radian')
866 case default
867 call mpp_error(fatal, "get_axis_latlon_data : Units for '// &
868 & trim(name)//' not recognised in file "//trim(fileobj%path))
869 end select
870
871 end subroutine get_axis_latlon_data_
872
873
874subroutine get_axis_level_data_(fileobj, name, level_data, level_type, vertical_indices)
875 type(FmsNetcdfFile_t), intent(in) :: fileobj
876 character(len=*), intent(in) :: name
877 real(FMS_INTP_KIND_), dimension(:), intent(out) :: level_data
878 integer, intent(out) :: level_type, vertical_indices
879 real(FMS_INTP_KIND_), dimension(:), allocatable :: alpha
880 integer :: n, nlev
881
882 integer, parameter :: lkind=fms_intp_kind_
883
884 if(variable_exists(fileobj, name)) then
885 call fms2_io_read_data(fileobj, name, level_data)
886 else
887 call mpp_error(fatal,'get_axis_level_data: variable '// &
888 & trim(name)//' does not exist in file '//trim(fileobj%path) )
889 endif
890 call get_variable_units(fileobj, name, units)
891 level_type = pressure
892 ! Convert to Pa
893 if( trim(adjustl(lowercase(chomp(units)))) == "mb" .or. trim(adjustl(lowercase(chomp(units)))) == "hpa") then
894 level_data = level_data * 100._lkind
895 endif
896 nlev = size(level_data(:))
897 sense = get_variable_sense(fileobj, name)
898! define the direction of the vertical data axis
899! switch index order if necessary so that indx 1 is at lowest pressure,
900! index nlev at highest pressure.
901 if( sense == 1 ) then
902 vertical_indices = increasing_upward
903 allocate (alpha(nlev))
904 do n = 1, nlev
905 alpha(n) = level_data(nlev-n+1)
906 end do
907 do n = 1, nlev
908 level_data(n) = alpha(n)
909 end do
910 deallocate (alpha)
911 else
912 vertical_indices = increasing_downward
913 endif
914
915
916end subroutine get_axis_level_data_
917
918
919!
920!---------------------------------------------------------------------
921!
922!> @brief cell_center2 receives the variables q1, q2, q3, and q4
923!! as inputs and returns e2.
924!!
925!! @param [in] <q1> No description
926!! @param [in] <q2> No description
927!! @param [in] <q3> No description
928!! @param [in] <q4> No description
929!! @param [out] <e2> No description
930 subroutine cell_center2_(q1, q2, q3, q4, e2)
931 real(FMS_INTP_KIND_) , intent(in ) :: q1(2), q2(2), q3(2), q4(2)
932 real(FMS_INTP_KIND_) , intent(out) :: e2(2)
933! Local
934 real(FMS_INTP_KIND_) :: p1(3), p2(3), p3(3), p4(3)
935 real(FMS_INTP_KIND_) :: ec(3,1)
936 real(FMS_INTP_KIND_) :: dd
937 integer k
938
939 call latlon2xyz(q1, p1)
940 call latlon2xyz(q2, p2)
941 call latlon2xyz(q3, p3)
942 call latlon2xyz(q4, p4)
943
944 do k=1,3
945 ec(k,1) = p1(k) + p2(k) + p3(k) + p4(k)
946 enddo
947 dd = sqrt( ec(1,1)**2 + ec(2,1)**2 + ec(3,1)**2 )
948
949 do k=1,3
950 ec(k,1) = ec(k,1) / dd
951 enddo
952
953 call cart_to_latlon(1, ec, e2(1:1), e2(2:2))
954
955 end subroutine cell_center2_
956!
957!---------------------------------------------------------------------
958!> @brief car_to_latlon receives the variables np, q, xs, and ys
959!! as inputs and returns q, xs, and ys.
960!!
961!! @param [in] <np> No description
962!! @param [inout] <q> No description
963!! @param [inout] <xs> No description
964!! @param [inout] <ys> No description
965 subroutine cart_to_latlon_(np, q, xs, ys)
966! vector version of cart_to_latlon1
967 integer, intent(in):: np
968 real(FMS_INTP_KIND_), intent(inout):: q(3,np)
969 real(FMS_INTP_KIND_), intent(inout):: xs(np), ys(np)
970! local
971 real(f_p), parameter:: esl=1.e-10_f_p
972 real(f_p):: p(3)
973 real(f_p):: dist, lat, lon
974 integer i,k
975
976 do i=1,np
977 do k=1,3
978 p(k) = real(q(k,i), f_p)
979 enddo
980 dist = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
981 do k=1,3
982 p(k) = p(k) / dist
983 enddo
984
985 if ( (abs(p(1))+abs(p(2))) < esl ) then
986 lon = 0._f_p
987 else
988 lon = atan2( p(2), p(1) ) ! range [-pi,pi]
989 endif
990
991 if ( lon < 0._f_p ) lon = 2._f_p*real(pi,f_p) + lon
992 lat = asin(p(3))
993
994 xs(i) = real(lon,fms_intp_kind_)
995 ys(i) = real(lat,fms_intp_kind_)
996! q Normalized:
997 do k=1,3
998 q(k,i) = real(p(k),fms_intp_kind_)
999 enddo
1000 enddo
1001
1002 end subroutine cart_to_latlon_
1003!
1004!---------------------------------------------------------------------
1005!> @brief latlon2xyz receives the variable p as input and returns e
1006!! as output in order to map (lon, lat) to (x,y,z).
1007!!
1008!! @param [in] <p> No description
1009!! @param [out] <e> No description
1010 subroutine latlon2xyz_(p, e)
1011!
1012! Routine to map (lon, lat) to (x,y,z)
1013!
1014 real(FMS_INTP_KIND_), intent(in) :: p(2)
1015 real(FMS_INTP_KIND_), intent(out):: e(3)
1016
1017 integer :: n
1018 real (f_p):: q(2)
1019 real (f_p):: e1, e2, e3
1020
1021 do n=1,2
1022 q(n) = real(p(n),f_p)
1023 enddo
1024
1025 e1 = cos(q(2)) * cos(q(1))
1026 e2 = cos(q(2)) * sin(q(1))
1027 e3 = sin(q(2))
1028!-----------------------------------
1029! Truncate to the desired precision:
1030!-----------------------------------
1031 e(1) = real(e1,fms_intp_kind_)
1032 e(2) = real(e2,fms_intp_kind_)
1033 e(3) = real(e3,fms_intp_kind_)
1034
1035 end subroutine latlon2xyz_
1036
1037!
1038!#######################################################################
1039!
1040!---------------------------------------------------------------------
1041!> @brief init_clim_diag is a routine to register diagnostic fields
1042!! for the climatology file. This routine calculates the domain
1043!! decompostion of the climatology fields for later export
1044!! through send_data. The ids created here are for column
1045!! burdens that will diagnose the vertical interpolation
1046!! routine.
1047!!
1048!! @param [inout] <clim_type> The interpolate type containing the
1049!! names of the fields in the climatology file
1050!! @param [in] <mod_axes> The axes of the model
1051!! @param [in] <init_time> The model initialization time
1052!!
1053!! @throw FATAL, "init_clim_diag : You must call interpolator_init before calling init_clim_diag"
1054!! @throw FATAL, "init_clim_diag : Trying to set up too many diagnostic fields for the climatology data"
1055subroutine init_clim_diag_(clim_type, mod_axes, init_time)
1056!
1057! Routine to register diagnostic fields for the climatology file.
1058! This routine calculates the domain decompostion of the climatology fields
1059! for later export through send_data.
1060! The ids created here are for column burdens that will diagnose the vertical interpolation routine.
1061! climo_diag_id : 'module_name = climo' is intended for use with the model vertical resolution.
1062! hinterp_id : 'module_name = 'hinterp' is intended for use with the climatology vertical resolution.
1063
1064! INTENT INOUT :
1065! clim_type : The interpolate type containing the names of the fields in the climatology file.
1066!
1067! INTENT IN :
1068! mod_axes : The axes of the model.
1069! init_time : The model initialization time.
1070!
1071type(interpolate_type), intent(inout) :: clim_type
1072integer , intent(in) :: mod_axes(:)
1073type(time_type) , intent(in) :: init_time
1074
1075integer :: axes(2),nxd,nyd,ndivs,i
1076type(domain2d) :: domain
1077integer :: domain_layout(2), iscomp, iecomp,jscomp,jecomp
1078
1079
1080if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) &
1081 call mpp_error(fatal, "init_clim_diag : You must call interpolator_init before calling init_clim_diag")
1082
1083
1084ndivs = mpp_npes()
1085nxd = size(clim_type%FMS_INTP_TYPE_%lon(:))
1086nyd = size(clim_type%FMS_INTP_TYPE_%lat(:))
1087
1088! Define the domain decomposition of the climatology file. This may be (probably is) different from the model domain.
1089call mpp_define_layout ((/1,nxd,1,nyd/), ndivs, domain_layout)
1090call mpp_define_domains((/1,nxd,1,nyd/),domain_layout, domain,xhalo=0,yhalo=0)
1091call mpp_get_compute_domain(domain, iscomp, iecomp, jscomp, jecomp)
1092axes(1) = diag_axis_init(clim_type%file_name(1:5)//'x',clim_type%FMS_INTP_TYPE_%lon,units='degrees',&
1093 cart_name='x',domain2=domain)
1094axes(2) = diag_axis_init(clim_type%file_name(1:5)//'y',clim_type%FMS_INTP_TYPE_%lat,units='degrees',&
1095 cart_name='y',domain2=domain)
1096clim_type%is = iscomp
1097clim_type%ie = iecomp
1098clim_type%js = jscomp
1099clim_type%je = jecomp
1100
1101!init_time = set_date(1980,1,1,0,0,0)
1102
1103if ((num_clim_diag + size(clim_type%field_name(:))) .gt. max_diag_fields ) &
1104 call mpp_error(fatal, "init_clim_diag : Trying to set up too many diagnostic fields for the climatology data")
1105do i=1,size(clim_type%field_name(:))
1106climo_diag_name(i+num_clim_diag) = clim_type%field_name(i)
1107climo_diag_id(i+num_clim_diag) = register_diag_field('climo',clim_type%field_name(i),axes(1:2),init_time,&
1108 'climo_'//clim_type%field_name(i), 'kg/kg', missing_value)
1109hinterp_id(i+num_clim_diag) = register_diag_field('hinterp',clim_type%field_name(i),mod_axes(1:2),init_time,&
1110 'interp_'//clim_type%field_name(i),'kg/kg' , missing_value)
1111enddo
1112! Total number of climatology diagnostics (num_clim_diag). This can be from multiple climatology
1113! fields with different spatial axes.
1114! It is simply a holder for the diagnostic indices.
1115num_clim_diag = num_clim_diag+size(clim_type%field_name(:))
1116
1117clim_diag_initialized = .true.
1118
1119end subroutine init_clim_diag_
1120
1121
1122
1123!
1124!---------------------------------------------------------------------
1125!> @brief obtain_interpolator_time_slices makes sure that the
1126!! appropriate time slices are available for interpolation on
1127!! this time step.
1128!!
1129!! @param [inout] <clim_type> The interpolate type previously defined
1130!! by a call to interpolator_init
1131!! @param [in] <Time> The model time that you wish to interpolate to
1132!!
1133!! @throw FATAL "interpolator_timeslice 1: file="
1134!! @throw FATAL "interpolator_timeslice 2: file="
1135!! @throw FATAL "interpolator_timeslice 3: file="
1136!! @throw FATAL "interpolator_timeslice 4: file="
1137!! @throw FATAL "interpolator_timeslice 5: file="
1138!! @throw FATAL "interpolator_timeslice : No data from the previous
1139!! climatology time but we have the next time. How did
1140!! this happen?"
1141subroutine obtain_interpolator_time_slices_ (clim_type, Time)
1142
1143! Makes sure that appropriate time slices are available for interpolation
1144! on this time step
1145!
1146! INTENT INOUT
1147! clim_type : The interpolate type previously defined by a call to interpolator_init
1148!
1149! INTENT IN
1150! Time : The model time that you wish to interpolate to.
1151
1152type(interpolate_type), intent(inout) :: clim_type
1153type(time_type) , intent(in) :: Time
1154
1155integer :: taum, taup
1156integer :: modyear, modmonth, modday, modhour, modminute, modsecond
1157integer :: climyear, climmonth, climday, climhour, climminute, climsecond
1158integer :: year1, month1, day, hour, minute, second
1159integer :: climatology, m
1160type(time_type) :: t_prev, t_next
1161type(time_type), dimension(2) :: month
1162integer :: indexm, indexp, yearm, yearp
1163integer :: i, n
1164character(len=256) :: err_msg
1165
1166integer, parameter :: lkind=fms_intp_kind_
1167
1168 if (clim_type%climatological_year) then
1169 !++lwh
1170 if (size(clim_type%time_slice) > 1) then
1171 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, &
1172 taum, taup, modtime=year, err_msg=err_msg )
1173 if(trim(err_msg) /= '') then
1174 call mpp_error(fatal,'interpolator_timeslice 1: '// &
1175 & trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1176 endif
1177 else
1178 taum = 1
1179 taup = 1
1180 clim_type%FMS_INTP_TYPE_%tweight = 0._lkind
1181 end if
1182 !--lwh
1183 else
1184 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, taum, taup, err_msg=err_msg )
1185 if(trim(err_msg) /= '') then
1186 call mpp_error(fatal,'interpolator_timeslice 2: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1187 endif
1188 endif
1189
1190 if(clim_type%TIME_FLAG .eq. bilinear ) then
1191 ! Check if delta-time is greater than delta of first two climatology time-slices.
1192 if ( (time - clim_type%time_slice(taum) ) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) .or. &
1193 (clim_type%time_slice(taup) - time) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) ) then
1194 ! The difference between the model time and the last climatology time-slice previous to the model time.
1195 ! We need 2 time levels.
1196 clim_type%itaum=0
1197 clim_type%itaup=0
1198 ! Assume this is monthly data. So we need to get the data applicable to the model date but substitute
1199 ! the climatology year into the appropriate place.
1200
1201
1202 ! We need to get the previous months data for the climatology year before
1203 ! and after the model year.
1204 call get_date(time, modyear, modmonth, modday, modhour, modminute, modsecond)
1205 call get_date(clim_type%time_slice(taum), climyear, climmonth, climday, climhour, climminute, climsecond)
1206
1207 climatology = 1
1208 do m = 1, size(clim_type%clim_times(:,:),2)
1209 !Assume here that a climatology is for 1 year and consists of 12 months starting in January.
1210 call get_date(clim_type%clim_times(1,m), year1, month1, day, hour, minute, second)
1211 if (year1 == climyear) climatology = m
1212 enddo
1213 do m = 1,12
1214 !Find which month we are trying to look at and set clim_date[mp] to the dates spanning that.
1215 call get_date(clim_type%clim_times(m,climatology), year1, month1, day, hour, minute, second)
1216 if ( month1 == modmonth ) then
1217!RSHBUGFX if ( modday <= day ) then
1218 if ( modday < day ) then
1219 indexm = m-1 ; indexp = m
1220 else
1221 indexm = m ; indexp = m+1
1222 endif
1223 endif
1224
1225 enddo
1226 if ( indexm == 0 ) then
1227 indexm = 12
1228 yearm = modyear - 1
1229 else
1230 yearm = modyear
1231 endif
1232 call get_date(clim_type%time_slice(indexm+(climatology-1)*12), &
1233 climyear, climmonth, climday, climhour, climminute, climsecond)
1234 month(1) = set_date(yearm, indexm, climday, climhour, climminute, climsecond)
1235 if ( indexp == 13 ) then
1236 indexp = 1
1237 yearp = modyear + 1
1238 else
1239 yearp = modyear
1240 endif
1241 call get_date(clim_type%time_slice(indexp+(climatology-1)*12), &
1242 climyear, climmonth, climday, climhour, climminute, climsecond)
1243 month(2) = set_date(yearp, indexp, climday, climhour, climminute, climsecond)
1244
1245 call time_interp(time, month, clim_type%FMS_INTP_TYPE_%tweight3, taum, taup, err_msg=err_msg ) ! tweight3 is
1246 !! the time weight between the months.
1247 !> protect against post-perth time_interp behavior below
1248 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight3 = 1._lkind
1249 if(trim(err_msg) /= '') then
1250 call mpp_error(fatal,'interpolator_timeslice 3: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1251 endif
1252
1253 month(1) = clim_type%time_slice(indexm+(climatology-1)*12)
1254 month(2) = clim_type%time_slice(indexm+climatology*12)
1255 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
1256 t_prev = set_date(yearm, climmonth, climday, climhour, climminute, climsecond)
1257 call time_interp(t_prev, month, clim_type%FMS_INTP_TYPE_%tweight1, taum, taup, err_msg=err_msg ) !tweight1 is
1258 !! the time weight between the climatology years.
1259 !> protect against post-perth time_interp behavior below
1260 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight1 = 1._lkind
1261 if(trim(err_msg) /= '') then
1262 call mpp_error(fatal,'interpolator_timeslice 4: '// &
1263 & trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1264 endif
1265
1266 month(1) = clim_type%time_slice(indexp+(climatology-1)*12)
1267 month(2) = clim_type%time_slice(indexp+climatology*12)
1268 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
1269 t_next = set_date(yearp, climmonth, climday, climhour, climminute, climsecond)
1270 call time_interp(t_next, month, clim_type%FMS_INTP_TYPE_%tweight2, taum, taup, err_msg=err_msg ) !tweight1 is
1271 !! the time weight between the climatology years.
1272 !> protect against post-perth time_interp behavior below
1273 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight2 = 1._lkind
1274 if(trim(err_msg) /= '') then
1275 call mpp_error(fatal,'interpolator_timeslice 5: '// &
1276 & trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1277 endif
1278
1279 if (indexm == clim_type%indexm(1) .and. &
1280 indexp == clim_type%indexp(1) .and. &
1281 climatology == clim_type%climatology(1)) then
1282 else
1283 clim_type%indexm(:) = indexm
1284 clim_type%indexp(:) = indexp
1285 clim_type%climatology(:) = climatology
1286 do i=1, size(clim_type%field_name(:))
1287 call read_data(clim_type,clim_type%field_name(i), &
1288 clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), &
1289 clim_type%indexm(i)+(clim_type%climatology(i)-1)*12,i,time)
1290! Read the data for the next month in the previous climatology.
1291 call read_data(clim_type,clim_type%field_name(i), &
1292 clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), &
1293 clim_type%indexp(i)+(clim_type%climatology(i)-1)*12,i,time)
1294 call read_data(clim_type,clim_type%field_name(i), &
1295 clim_type%FMS_INTP_TYPE_%pmon_nyear(:,:,:,i), &
1296 clim_type%indexm(i)+clim_type%climatology(i)*12,i,time)
1297 call read_data(clim_type,clim_type%field_name(i), &
1298 clim_type%FMS_INTP_TYPE_%nmon_nyear(:,:,:,i), &
1299 clim_type%indexp(i)+clim_type%climatology(i)*12,i,time)
1300 end do
1301 endif
1302
1303
1304 else ! We are within a climatology data set
1305
1306 do i=1, size(clim_type%field_name(:))
1307 if (taum /= clim_type%time_init(i,1) .or. &
1308 taup /= clim_type%time_init(i,2) ) then
1309
1310
1311 call read_data(clim_type,clim_type%field_name(i), &
1312 clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), taum,i,time)
1313! Read the data for the next month in the previous climatology.
1314 call read_data(clim_type,clim_type%field_name(i), &
1315 clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), taup,i,time)
1316 clim_type%time_init(i,1) = taum
1317 clim_type%time_init(i,2) = taup
1318 endif
1319 end do
1320! clim_type%pmon_nyear = 0.0
1321! clim_type%nmon_nyear = 0.0
1322
1323! set to zero so when next return to bilinear section will be sure to
1324! have proper data (relevant when running fixed_year case for more than
1325! one year in a single job)
1326 clim_type%indexm(:) = 0
1327 clim_type%indexp(:) = 0
1328 clim_type%climatology(:) = 0
1329
1330
1331! clim_type%tweight3 = 0.0 ! This makes [pn]mon_nyear irrelevant. Set them to 0 to test.
1332 clim_type%FMS_INTP_TYPE_%tweight1 = 0.0_lkind
1333 clim_type%FMS_INTP_TYPE_%tweight2 = 0.0_lkind
1334 clim_type%FMS_INTP_TYPE_%tweight3 = clim_type%FMS_INTP_TYPE_%tweight
1335 endif
1336 endif !(BILINEAR)
1337
1338 if(clim_type%TIME_FLAG .eq. linear .and. &
1339 (.not. read_all_on_init) ) then
1340! We need 2 time levels. Check we have the correct data.
1341 clim_type%itaum=0
1342 clim_type%itaup=0
1343 do n=1,size(clim_type%time_init,2)
1344 if (clim_type%time_init(1,n) .eq. taum ) clim_type%itaum = n
1345 if (clim_type%time_init(1,n) .eq. taup ) clim_type%itaup = n
1346 enddo
1347
1348 if (clim_type%itaum.eq.0 .and. clim_type%itaup.eq.0) then
1349!Neither time is set so we need to read 2 time slices.
1350!Set up
1351! field(:,:,:,1) as the previous time slice.
1352! field(:,:,:,2) as the next time slice.
1353 do i=1, size(clim_type%field_name(:))
1354 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,1,i), taum,i,time)
1355 clim_type%time_init(i,1) = taum
1356 clim_type%itaum = 1
1357 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,2,i), taup,i,time)
1358 clim_type%time_init(i,2) = taup
1359 clim_type%itaup = 2
1360 end do
1361 endif ! clim_type%itaum.eq.clim_type%itaup.eq.0
1362 if (clim_type%itaum.eq.0 .and. clim_type%itaup.ne.0) then
1363! Can't think of a situation where we would have the next time level but not the previous.
1364 call mpp_error(fatal,'interpolator_timeslice : No data from the previous climatology time &
1365 & but we have the next time. How did this happen?')
1366 endif
1367 if (clim_type%itaum.ne.0 .and. clim_type%itaup.eq.0) then
1368!We have the previous time step but not the next time step data
1369 clim_type%itaup = 1
1370 if (clim_type%itaum .eq. 1 ) clim_type%itaup = 2
1371 do i=1, size(clim_type%field_name(:))
1372 call read_data(clim_type,clim_type%field_name(i), &
1373 clim_type%FMS_INTP_TYPE_%data5d(:,:,:,clim_type%itaup,i), taup,i, time)
1374 clim_type%time_init(i,clim_type%itaup)=taup
1375 end do
1376 endif
1377
1378
1379 endif! TIME_FLAG
1380
1381 clim_type%separate_time_vary_calc = .true.
1382
1383!-------------------------------------------------------------------
1384
1385
1386 end subroutine obtain_interpolator_time_slices_
1387
1388
1389!#####################################################################
1390!
1391!---------------------------------------------------------------------
1392!> @brief interpolator_4D receives a field name as input and
1393!! interpolates the field to model a 4D grid and time axis.
1394!!
1395!! @param [inout] <clim_type> The interpolate type previously defined by a call to interpolator_init
1396!! @param [in] <field_name> The name of a field that you wish to interpolate
1397!! @param [in] <Time> The model time that you wish to interpolate to
1398!! @param [in] <phalf> The half level model pressure field
1399!! @param [in] <is> OPTIONAL: Index for the physics window
1400!! @param [in] <js> OPTIONAL: Index for the physics window
1401!! @param [out] <interp_data> The model fields with the interpolated climatology data
1402!! @param [out] <clim_units> OPTIONAL: The units of field_name
1403!!
1404!! @throw FATAL "interpolator_4D : You must call interpolator_init
1405!! before calling interpolator"
1406!! @throw FATAL "interpolator_mod: cannot use 4D interface to interpolator for this file"
1407!! @throw FATAL "interpolator_4D 1: file="
1408!! @throw FATAL "interpolator_4D 2: file="
1409!! @throw FATAL "interpolator_4D 3: file="
1410!! @throw FATAL "interpolator_4D 4: file="
1411!! @throw FATAL "interpolator_4D 5: file="
1412!! @throw FATAL "interpolator_3D : No data from the previous climatology
1413!! time but we have the next time. How did this happen?"
1414!! @throw NOTE "Interpolator: model surface pressure is greater than
1415!! climatology surface pressure for "
1416!! @throw NOTE "Interpolator: model top pressure is less than
1417!! climatology top pressure for "
1418!! @throw FATAL "Interpolator: the field name is not contained in this
1419!! intepolate_type: "
1420subroutine interpolator_4d_(clim_type, Time, phalf, interp_data, &
1421 field_name, is,js, clim_units)
1422!
1423! Return 4-D field interpolated to model grid and time
1424!
1425! INTENT INOUT
1426! clim_type : The interpolate type previously defined by a call to interpolator_init
1427!
1428! INTENT IN
1429! field_name : The name of a field that you wish to interpolate.
1430! all variables within this interpolate_type variable
1431! will be interpolated on this call. field_name may
1432! be any one of the variables.
1433! Time : The model time that you wish to interpolate to.
1434! phalf : The half level model pressure field.
1435! is, js : The indices of the physics window.
1436!
1437! INTENT OUT
1438! interp_data : The model fields with the interpolated climatology data.
1439! clim_units : The units of field_name
1440!
1441type(interpolate_type), intent(inout) :: clim_type
1442character(len=*) , intent(in) :: field_name
1443type(time_type) , intent(in) :: Time
1444real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf
1445real(FMS_INTP_KIND_), dimension(:,:,:,:), intent(out) :: interp_data
1446integer , intent(in) , optional :: is,js
1447character(len=*) , intent(out), optional :: clim_units
1448integer :: taum, taup, ilon
1449real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),size(interp_data,2),&
1450 size(clim_type%FMS_INTP_TYPE_%levs(:)),size(clim_type%field_name(:)))
1451real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2))
1452real(FMS_INTP_KIND_) :: col_data(size(interp_data,1),size(interp_data,2), &
1453 size(clim_type%field_name(:)))
1454real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:)))
1455integer :: istart,iend,jstart,jend
1456logical :: result, found
1457logical :: found_field=.false.
1458integer :: modyear, modmonth, modday, modhour, modminute, modsecond
1459integer :: climyear, climmonth, climday, climhour, climminute, climsecond
1460integer :: year1, month1, day, hour, minute, second
1461integer :: climatology, m
1462type(time_type) :: t_prev, t_next
1463type(time_type), dimension(2) :: month
1464integer :: indexm, indexp, yearm, yearp
1465integer :: i, j, k, n
1466character(len=256) :: err_msg
1467
1468integer, parameter :: lkind=fms_intp_kind_
1469
1470if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) &
1471 call mpp_error(fatal, "interpolator_4D : You must call interpolator_init before calling interpolator")
1472
1473 do n=2,size(clim_type%field_name(:))
1474 if (clim_type%vert_interp(n) /= clim_type%vert_interp(n-1) .or. &
1475 clim_type%out_of_bounds(n) /= clim_type%out_of_bounds(n-1)) then
1476 if (mpp_pe() == mpp_root_pe() ) then
1477 print *, 'processing file ' // trim(clim_type%file_name)
1478 endif
1479 call mpp_error (fatal, 'interpolator_mod: &
1480 &cannot use 4D interface to interpolator for this file')
1481 endif
1482 end do
1483
1484
1485
1486
1487istart = 1
1488if (present(is)) istart = is
1489iend = istart - 1 + size(interp_data,1)
1490
1491jstart = 1
1492if (present(js)) jstart = js
1493jend = jstart - 1 + size(interp_data,2)
1494
1495 do i= 1,size(clim_type%field_name(:))
1496!!++lwh
1497 if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then
1498!--lwh
1499 found_field=.true.
1500 exit
1501 endif
1502end do
1503 i = 1
1504
1505 if(present(clim_units)) then
1506 call get_variable_units(clim_type%fileobj, clim_type%field_name(i), clim_units)
1507 endif
1508
1509
1510!----------------------------------------------------------------------
1511! skip the time interpolation portion of this routine if subroutine
1512! obtain_interpolator_time_slices has already been called on this
1513! stewp for this interpolate_type variable.
1514!----------------------------------------------------------------------
1515
1516if ( .not. clim_type%separate_time_vary_calc) then
1517! print *, 'TIME INTERPOLATION NOT SEPARATED 4d--', &
1518! trim(clim_type%file_name), mpp_pe()
1519
1520 if (clim_type%climatological_year) then
1521!++lwh
1522 if (size(clim_type%time_slice) > 1) then
1523 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, &
1524 taum, taup, modtime=year, err_msg=err_msg )
1525 if(trim(err_msg) /= '') then
1526 call mpp_error(fatal,'interpolator_4D 1: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1527 endif
1528 else
1529 taum = 1
1530 taup = 1
1531 clim_type%FMS_INTP_TYPE_%tweight = 0._lkind
1532 end if
1533!--lwh
1534 else
1535 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, taum, taup, err_msg=err_msg )
1536 if(trim(err_msg) /= '') then
1537 call mpp_error(fatal,'interpolator_4D 2: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1538 endif
1539 endif
1540
1541 if(clim_type%TIME_FLAG .eq. bilinear ) then
1542 ! Check if delta-time is greater than delta of first two climatology time-slices.
1543 if ( (time - clim_type%time_slice(taum) ) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) .or. &
1544 (clim_type%time_slice(taup) - time) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) ) then
1545 ! The difference between the model time and the last climatology time-slice previous to the model time.
1546 ! We need 2 time levels.
1547 clim_type%itaum=0
1548 clim_type%itaup=0
1549 ! Assume this is monthly data. So we need to get the data applicable to the model date but substitute
1550 ! the climatology year into the appropriate place.
1551
1552
1553 ! We need to get the previous months data for the climatology year before
1554 ! and after the model year.
1555 call get_date(time, modyear, modmonth, modday, modhour, modminute, modsecond)
1556 call get_date(clim_type%time_slice(taum), climyear, climmonth, climday, climhour, climminute, climsecond)
1557
1558 climatology = 1
1559 do m = 1, size(clim_type%clim_times(:,:),2)
1560 !Assume here that a climatology is for 1 year and consists of 12 months starting in January.
1561 call get_date(clim_type%clim_times(1,m), year1, month1, day, hour, minute, second)
1562 if (year1 == climyear) climatology = m
1563 enddo
1564 do m = 1,12
1565 !Find which month we are trying to look at and set clim_date[mp] to the dates spanning that.
1566 call get_date(clim_type%clim_times(m,climatology), year1, month1, day, hour, minute, second)
1567 if ( month1 == modmonth ) then
1568!RSHBUGFX if ( modday <= day ) then
1569 if ( modday < day ) then
1570 indexm = m-1 ; indexp = m
1571 else
1572 indexm = m ; indexp = m+1
1573 endif
1574 endif
1575
1576 enddo
1577 if ( indexm == 0 ) then
1578 indexm = 12
1579 yearm = modyear - 1
1580 else
1581 yearm = modyear
1582 endif
1583 call get_date(clim_type%time_slice(indexm+(climatology-1)*12), &
1584 climyear, climmonth, climday, climhour, climminute, climsecond)
1585 month(1) = set_date(yearm, indexm, climday, climhour, climminute, climsecond)
1586 if ( indexp == 13 ) then
1587 indexp = 1
1588 yearp = modyear + 1
1589 else
1590 yearp = modyear
1591 endif
1592 call get_date(clim_type%time_slice(indexp+(climatology-1)*12), &
1593 climyear, climmonth, climday, climhour, climminute, climsecond)
1594 month(2) = set_date(yearp, indexp, climday, climhour, climminute, climsecond)
1595
1596 call time_interp(time, month, clim_type%FMS_INTP_TYPE_%tweight3, taum, taup, err_msg=err_msg ) ! tweight3 is
1597 !! the time weight between the months.
1598 !> protect against post-perth time_interp behavior below
1599 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight3 = 1._lkind
1600 if(trim(err_msg) /= '') then
1601 call mpp_error(fatal,'interpolator_4D 3: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1602 endif
1603
1604 month(1) = clim_type%time_slice(indexm+(climatology-1)*12)
1605 month(2) = clim_type%time_slice(indexm+climatology*12)
1606 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
1607 t_prev = set_date(yearm, climmonth, climday, climhour, climminute, climsecond)
1608 call time_interp(t_prev, month, clim_type%FMS_INTP_TYPE_%tweight1, taum, taup, err_msg=err_msg ) !tweight1 is
1609 !! the time weight between the climatology years.
1610 !> protect against post-perth time_interp behavior below
1611 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight1 = 1._lkind
1612 if(trim(err_msg) /= '') then
1613 call mpp_error(fatal,'interpolator_4D 4: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1614 endif
1615 month(1) = clim_type%time_slice(indexp+(climatology-1)*12)
1616 month(2) = clim_type%time_slice(indexp+climatology*12)
1617 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
1618 t_next = set_date(yearp, climmonth, climday, climhour, climminute, climsecond)
1619 call time_interp(t_next, month, clim_type%FMS_INTP_TYPE_%tweight2, taum, taup, err_msg=err_msg ) !tweight1 is
1620 !! the time weight between the climatology years.
1621 !> protect against post-perth time_interp behavior below
1622 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight2 = 1._lkind
1623 if(trim(err_msg) /= '') then
1624 call mpp_error(fatal,'interpolator_4D 5: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1625 endif
1626
1627 if (indexm == clim_type%indexm(1) .and. &
1628 indexp == clim_type%indexp(1) .and. &
1629 climatology == clim_type%climatology(1)) then
1630 else
1631 clim_type%indexm(:) = indexm
1632 clim_type%indexp(:) = indexp
1633 clim_type%climatology(:) = climatology
1634 do i=1, size(clim_type%field_name(:))
1635 call read_data(clim_type,clim_type%field_name(i), &
1636 clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), &
1637 clim_type%indexm(i)+(clim_type%climatology(i)-1)*12,i,time)
1638! Read the data for the next month in the previous climatology.
1639 call read_data(clim_type,clim_type%field_name(i), &
1640 clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), &
1641 clim_type%indexp(i)+(clim_type%climatology(i)-1)*12,i,time)
1642 call read_data(clim_type,clim_type%field_name(i), &
1643 clim_type%FMS_INTP_TYPE_%pmon_nyear(:,:,:,i), &
1644 clim_type%indexm(i)+clim_type%climatology(i)*12,i,time)
1645 call read_data(clim_type,clim_type%field_name(i), &
1646 clim_type%FMS_INTP_TYPE_%nmon_nyear(:,:,:,i), &
1647 clim_type%indexp(i)+clim_type%climatology(i)*12,i,time)
1648 end do
1649 endif
1650
1651
1652
1653 else ! We are within a climatology data set
1654
1655 do i=1, size(clim_type%field_name(:))
1656 if (taum /= clim_type%time_init(i,1) .or. &
1657 taup /= clim_type%time_init(i,2) ) then
1658
1659
1660 call read_data(clim_type,clim_type%field_name(i), &
1661 clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), taum,i,time)
1662! Read the data for the next month in the previous climatology.
1663 call read_data(clim_type,clim_type%field_name(i), &
1664 clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), taup,i,time)
1665 clim_type%time_init(i,1) = taum
1666 clim_type%time_init(i,2) = taup
1667 endif
1668 end do
1669! clim_type%pmon_nyear = 0.0
1670! clim_type%nmon_nyear = 0.0
1671
1672! set to zero so when next return to bilinear section will be sure to
1673! have proper data (relevant when running fixed_year case for more than
1674! one year in a single job)
1675 clim_type%indexm(:) = 0
1676 clim_type%indexp(:) = 0
1677 clim_type%climatology(:) = 0
1678
1679
1680! clim_type%tweight3 = 0.0 ! This makes [pn]mon_nyear irrelevant. Set them to 0 to test.
1681 clim_type%FMS_INTP_TYPE_%tweight1 = 0.0_lkind
1682 clim_type%FMS_INTP_TYPE_%tweight2 = 0.0_lkind
1683 clim_type%FMS_INTP_TYPE_%tweight3 = clim_type%FMS_INTP_TYPE_%tweight
1684 endif
1685 endif !(BILINEAR)
1686
1687 if(clim_type%TIME_FLAG .eq. linear .and. &
1688 (.not. read_all_on_init) ) then
1689! We need 2 time levels. Check we have the correct data.
1690 clim_type%itaum=0
1691 clim_type%itaup=0
1692 do n=1,size(clim_type%time_init,2)
1693 if (clim_type%time_init(1,n) .eq. taum ) clim_type%itaum = n
1694 if (clim_type%time_init(1,n) .eq. taup ) clim_type%itaup = n
1695 enddo
1696
1697 if (clim_type%itaum.eq.0 .and. clim_type%itaup.eq.0) then
1698!Neither time is set so we need to read 2 time slices.
1699!Set up
1700! field(:,:,:,1) as the previous time slice.
1701! field(:,:,:,2) as the next time slice.
1702 do i=1, size(clim_type%field_name(:))
1703 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,1,i), taum,i,time)
1704 clim_type%time_init(i,1) = taum
1705 clim_type%itaum = 1
1706 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,2,i), taup,i,time)
1707 clim_type%time_init(i,2) = taup
1708 clim_type%itaup = 2
1709 end do
1710 endif ! clim_type%itaum.eq.clim_type%itaup.eq.0
1711 if (clim_type%itaum.eq.0 .and. clim_type%itaup.ne.0) then
1712! Can't think of a situation where we would have the next time level but not the previous.
1713 call mpp_error(fatal,'interpolator_3D : No data from the previous climatology time &
1714 & but we have the next time. How did this happen?')
1715 endif
1716 if (clim_type%itaum.ne.0 .and. clim_type%itaup.eq.0) then
1717!We have the previous time step but not the next time step data
1718 clim_type%itaup = 1
1719 if (clim_type%itaum .eq. 1 ) clim_type%itaup = 2
1720 do i=1, size(clim_type%field_name(:))
1721 call read_data(clim_type,clim_type%field_name(i), &
1722 clim_type%FMS_INTP_TYPE_%data5d(:,:,:,clim_type%itaup,i), taup,i, time)
1723 clim_type%time_init(i,clim_type%itaup)=taup
1724 end do
1725 endif
1726
1727
1728 endif! TIME_FLAG
1729
1730
1731endif ! (.not. separate_time_vary_calc)
1732
1733
1734select case(clim_type%TIME_FLAG)
1735 case (linear)
1736 do n=1, size(clim_type%field_name(:))
1737 hinterp_data(:,:,:,n) = (1._lkind-clim_type%FMS_INTP_TYPE_%tweight)* &
1738 clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,clim_type%itaum,n) + &
1739 clim_type%FMS_INTP_TYPE_%tweight* &
1740 clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,clim_type%itaup,n)
1741 end do
1742! case (SEASONAL)
1743! Do sine fit to data at this point
1744 case (bilinear)
1745 do n=1, size(clim_type%field_name(:))
1746 hinterp_data(:,:,:,n) = (1._lkind-clim_type%FMS_INTP_TYPE_%tweight1)*&
1747 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight3)* &
1748 clim_type%FMS_INTP_TYPE_%pmon_pyear(istart:iend,jstart:jend,:,n) + &
1749 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight2)*clim_type%FMS_INTP_TYPE_%tweight3* &
1750 clim_type%FMS_INTP_TYPE_%nmon_pyear(istart:iend,jstart:jend,:,n) + &
1751 clim_type%FMS_INTP_TYPE_%tweight1* (1._lkind-clim_type%FMS_INTP_TYPE_%tweight3)*&
1752 clim_type%FMS_INTP_TYPE_%pmon_nyear(istart:iend,jstart:jend,:,n) + &
1753 clim_type%FMS_INTP_TYPE_%tweight2* clim_type%FMS_INTP_TYPE_%tweight3* &
1754 clim_type%FMS_INTP_TYPE_%nmon_nyear(istart:iend,jstart:jend,:,n)
1755
1756 end do
1757
1758end select
1759
1760select case(clim_type%level_type)
1761 case(pressure)
1762 p_fact = 1.0_lkind
1763 case(sigma)
1764 p_fact = maxval(phalf,3)! max pressure in the column !(:,:,size(phalf,3))
1765end select
1766
1767col_data(:,:,:)=0.0_lkind
1768 do i= 1, size(clim_type%field_name(:))
1769
1770select case(clim_type%mr(i))
1771 case(no_conv)
1772 do k = 1,size(hinterp_data,3)
1773 col_data(:,:,i) = col_data(:,:,i) + hinterp_data(:,:,k,i)* &
1774 (clim_type%FMS_INTP_TYPE_%halflevs(k+1)-clim_type%FMS_INTP_TYPE_%halflevs(k))/real(grav,fms_intp_kind_)
1775 enddo
1776
1777 case(kg_m2)
1778 do k = 1,size(hinterp_data,3)
1779 col_data(:,:,i) = col_data(:,:,i) + hinterp_data(:,:,k,i)
1780 hinterp_data(:,:,k,i) = hinterp_data(:,:,k,i)/ &
1781 ((clim_type%FMS_INTP_TYPE_%halflevs(k+1)-clim_type%FMS_INTP_TYPE_%halflevs(k))*p_fact)
1782 enddo
1783end select
1784 enddo
1785
1786 do i= 1, size(clim_type%field_name(:))
1787found = .false.
1788do j = 1,size(climo_diag_name(:))
1789 if ( trim(adjustl(lowercase(climo_diag_name(j)))) .eq. trim(adjustl(lowercase(clim_type%field_name(i))))) then
1790 found = .true.
1791 exit
1792 endif
1793enddo
1794
1795if (found) then
1796 if (hinterp_id(j) > 0 ) then
1797 result = send_data(hinterp_id(j),col_data(:,:,i),time,is_in=istart,js_in=jstart)
1798 endif
1799endif
1800
1801 end do
1802
1803 i = 1
1804
1805!++lwh
1806do j = 1, size(phalf,2)
1807 do ilon=1,size(phalf,1)
1808 pclim = p_fact(ilon,j)*clim_type%FMS_INTP_TYPE_%halflevs
1809 if ( maxval(phalf(ilon,j,:)) > maxval(pclim) ) then
1810 if (verbose > 3) then
1811 call mpp_error(note,"Interpolator: model surface pressure&
1812 & is greater than climatology surface pressure for "&
1813 // trim(clim_type%file_name))
1814 endif
1815 select case(clim_type%out_of_bounds(i))
1816 case(constant)
1817 pclim( maxloc(pclim) ) = maxval( phalf(ilon,j,:) )
1818! case(ZERO)
1819! pclim( maxloc(pclim)) = 0
1820 end select
1821 endif
1822 if ( minval(phalf(ilon,j,:)) < minval(pclim) ) then
1823 if (verbose > 3) then
1824 call mpp_error(note,"Interpolator: model top pressure&
1825 & is less than climatology top pressure for "&
1826 & // trim(clim_type%file_name))
1827 endif
1828 select case(clim_type%out_of_bounds(i))
1829 case(constant)
1830 pclim( minloc(pclim) ) = minval( phalf(ilon,j,:) )
1831! case(ZERO)
1832! pclim( maxloc(pclim)) = 0
1833 end select
1834 endif
1835 select case(clim_type%vert_interp(i))
1836 case(interp_weighted_p)
1837 call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,:),interp_data(ilon,j,:,:))
1838 case(interp_linear_p)
1839 do n=1, size(clim_type%field_name(:))
1840 call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,n),interp_data(ilon,j,:,n))
1841 end do
1842! case(INTERP_LOG)
1843 end select
1844 enddo
1845enddo
1846
1847!--lwh
1848 do i= 1, size(clim_type%field_name(:))
1849
1850select case(clim_type%mr(i))
1851 case(kg_m2)
1852 do k = 1,size(interp_data,3)
1853 interp_data(:,:,k,i) = interp_data(:,:,k,i)*(phalf(:,:,k+1)-phalf(:,:,k))
1854 enddo
1855end select
1856
1857 end do
1858
1859if( .not. found_field) then !field name is not in interpolator file.ERROR.
1860 call mpp_error(fatal,"Interpolator: the field name is not contained in this &
1861 &intepolate_type: "//trim(field_name))
1862endif
1863end subroutine interpolator_4d_
1864!
1865!#######################################################################
1866!#######################################################################
1867!
1868!---------------------------------------------------------------------
1869!> @brief interpolator_3D receives a field name as input and
1870!! interpolates the field to model a 3D grid and time axis.
1871!!
1872!! @param [inout] <clim_type> The interpolate type previously defined by a call to interpolator_init
1873!! @param [in] <field_name> The name of a field that you wish to interpolate
1874!! @param [in] <Time> The model time that you wish to interpolate to
1875!! @param [in] <phalf> The half level model pressure field
1876!! @param [in] <is> OPTIONAL: Index for the physics window
1877!! @param [in] <js> OPTIONAL: Index for the physics window
1878!! @param [out] <interp_data> The model fields with the interpolated climatology data
1879!! @param [out] <clim_units> OPTIONAL: The units of field_name
1880!!
1881!! @throw FATAL "interpolator_3D : You must call interpolator_init
1882!! before calling interpolator"
1883!! @throw FATAL "interpolator_3D 1: file="
1884!! @throw FATAL "interpolator_3D 2: file="
1885!! @throw FATAL "interpolator_3D 3: file="
1886!! @throw FATAL "interpolator_3D 4: file="
1887!! @throw FATAL "interpolator_3D 5: file="
1888!! @throw FATAL "interpolator_3D : No data from the previous climatology
1889!! time but we have the next time. How did this happen?"
1890!! @throw NOTE "Interpolator: model surface pressure is greater than
1891!! climatology surface pressure for "
1892!! @throw NOTE "Interpolator: model top pressure is less than
1893!! climatology top pressure for "
1894!! @throw FATAL "Interpolator: the field name is not contained in this
1895!! intepolate_type: "
1896subroutine interpolator_3d_(clim_type, Time, phalf, interp_data,field_name, is,js, clim_units)
1897!
1898! Return 3-D field interpolated to model grid and time
1899!
1900! INTENT INOUT
1901! clim_type : The interpolate type previously defined by a call to interpolator_init
1902!
1903! INTENT IN
1904! field_name : The name of the field that you wish to interpolate.
1905! Time : The model time that you wish to interpolate to.
1906! phalf : The half level model pressure field.
1907! is, js : The indices of the physics window.
1908!
1909! INTENT OUT
1910! interp_data : The model field with the interpolated climatology data.
1911! clim_units : The units of field_name
1912!
1913type(interpolate_type), intent(inout) :: clim_type
1914character(len=*) , intent(in) :: field_name
1915type(time_type) , intent(in) :: Time
1916real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf
1917real(FMS_INTP_KIND_), dimension(:,:,:), intent(out) :: interp_data
1918integer , intent(in) , optional :: is,js
1919character(len=*) , intent(out), optional :: clim_units
1920integer :: taum, taup, ilon
1921real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),&
1922 size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:)))
1923real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2))
1924real(FMS_INTP_KIND_) :: col_data(size(interp_data,1),size(interp_data,2))
1925real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:)))
1926integer :: istart,iend,jstart,jend
1927logical :: result, found
1928logical :: found_field=.false.
1929integer :: modyear, modmonth, modday, modhour, modminute, modsecond
1930integer :: climyear, climmonth, climday, climhour, climminute, climsecond
1931integer :: year1, month1, day, hour, minute, second
1932integer :: climatology, m
1933type(time_type) :: t_prev, t_next
1934type(time_type), dimension(2) :: month
1935integer :: indexm, indexp, yearm, yearp
1936integer :: i, j, k, n
1937character(len=256) :: err_msg
1938
1939integer, parameter :: lkind=fms_intp_kind_
1940
1941if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) &
1942 call mpp_error(fatal, "interpolator_3D : You must call interpolator_init before calling interpolator")
1943
1944istart = 1
1945if (present(is)) istart = is
1946iend = istart - 1 + size(interp_data,1)
1947
1948jstart = 1
1949if (present(js)) jstart = js
1950jend = jstart - 1 + size(interp_data,2)
1951
1952do i= 1,size(clim_type%field_name(:))
1953!++lwh
1954 if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then
1955!--lwh
1956 found_field=.true.
1957 if(present(clim_units)) then
1958 call get_variable_units(clim_type%fileobj, clim_type%field_name(i), clim_units)
1959 clim_units = chomp(clim_units)
1960 endif
1961
1962!----------------------------------------------------------------------
1963! skip the time interpolation portion of this routine if subroutine
1964! obtain_interpolator_time_slices has already been called on this
1965! stewp for this interpolate_type variable.
1966!----------------------------------------------------------------------
1967
1968if ( .not. clim_type%separate_time_vary_calc) then
1969! print *, 'TIME INTERPOLATION NOT SEPARATED 3d--', &
1970! trim(clim_type%file_name), mpp_pe()
1971 if (clim_type%climatological_year) then
1972!++lwh
1973 if (size(clim_type%time_slice) > 1) then
1974 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, &
1975 taum, taup, modtime=year, err_msg=err_msg )
1976 if(trim(err_msg) /= '') then
1977 call mpp_error(fatal,'interpolator_3D 1: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1978 endif
1979 else
1980 taum = 1
1981 taup = 1
1982 clim_type%FMS_INTP_TYPE_%tweight = 0._lkind
1983 end if
1984!--lwh
1985 else
1986 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, taum, taup, err_msg=err_msg )
1987 if(trim(err_msg) /= '') then
1988 call mpp_error(fatal,'interpolator_3D 2: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
1989 endif
1990 endif
1991
1992! if(clim_type%TIME_FLAG .ne. LINEAR ) then
1993 if(clim_type%TIME_FLAG .ne. linear .or. read_all_on_init ) then
1994 clim_type%itaum=taum
1995 clim_type%itaup=taup
1996 endif
1997
1998 if(clim_type%TIME_FLAG .eq. bilinear ) then
1999 ! Check if delta-time is greater than delta of first two climatology time-slices.
2000 if ( (time - clim_type%time_slice(taum) ) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) .or. &
2001 (clim_type%time_slice(taup) - time) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) ) then
2002 ! The difference between the model time and the last climatology time-slice previous to the model time.
2003 ! We need 2 time levels.
2004 clim_type%itaum=0
2005 clim_type%itaup=0
2006 ! Assume this is monthly data. So we need to get the data applicable to the model date but substitute
2007 ! the climatology year into the appropriate place.
2008
2009
2010 ! We need to get the previous months data for the climatology year before
2011 ! and after the model year.
2012 call get_date(time, modyear, modmonth, modday, modhour, modminute, modsecond)
2013 call get_date(clim_type%time_slice(taum), climyear, climmonth, climday, climhour, climminute, climsecond)
2014
2015 climatology = 1
2016 do m = 1, size(clim_type%clim_times(:,:),2)
2017 !Assume here that a climatology is for 1 year and consists of 12 months starting in January.
2018 call get_date(clim_type%clim_times(1,m), year1, month1, day, hour, minute, second)
2019 if (year1 == climyear) climatology = m
2020 enddo
2021 do m = 1,12
2022 !Find which month we are trying to look at and set clim_date[mp] to the dates spanning that.
2023 call get_date(clim_type%clim_times(m,climatology), year1, month1, day, hour, minute, second)
2024 if ( month1 == modmonth ) then
2025!RSHBUGFX if ( modday <= day ) then
2026 if ( modday < day ) then
2027 indexm = m-1 ; indexp = m
2028 else
2029 indexm = m ; indexp = m+1
2030 endif
2031 endif
2032
2033 enddo
2034 if ( indexm == 0 ) then
2035 indexm = 12
2036 yearm = modyear - 1
2037 else
2038 yearm = modyear
2039 endif
2040 call get_date(clim_type%time_slice(indexm+(climatology-1)*12), &
2041 climyear, climmonth, climday, climhour, climminute, climsecond)
2042 month(1) = set_date(yearm, indexm, climday, climhour, climminute, climsecond)
2043 if ( indexp == 13 ) then
2044 indexp = 1
2045 yearp = modyear + 1
2046 else
2047 yearp = modyear
2048 endif
2049 call get_date(clim_type%time_slice(indexp+(climatology-1)*12), &
2050 climyear, climmonth, climday, climhour, climminute, climsecond)
2051 month(2) = set_date(yearp, indexp, climday, climhour, climminute, climsecond)
2052
2053 call time_interp(time, month, clim_type%FMS_INTP_TYPE_%tweight3, taum, taup, err_msg=err_msg ) ! tweight3 is
2054 !! the time weight between the months.
2055 if (taum==2 .and. taup==2) &
2056 clim_type%FMS_INTP_TYPE_%tweight3 = 1._lkind ! protect against post-perth time_interp behavior
2057 if(trim(err_msg) /= '') then
2058 call mpp_error(fatal,'interpolator_3D 3: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2059 endif
2060
2061 month(1) = clim_type%time_slice(indexm+(climatology-1)*12)
2062 month(2) = clim_type%time_slice(indexm+climatology*12)
2063 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
2064 t_prev = set_date(yearm, climmonth, climday, climhour, climminute, climsecond)
2065 call time_interp(t_prev, month, clim_type%FMS_INTP_TYPE_%tweight1, taum, taup, err_msg=err_msg ) !tweight1 is
2066 !! the time weight between the climatology years.
2067 if (taum==2 .and. taup==2) &
2068 clim_type%FMS_INTP_TYPE_%tweight1 = 1._lkind ! protect against post-perth time_interp behavior
2069 if(trim(err_msg) /= '') then
2070 call mpp_error(fatal,'interpolator_3D 4: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2071 endif
2072
2073 month(1) = clim_type%time_slice(indexp+(climatology-1)*12)
2074 month(2) = clim_type%time_slice(indexp+climatology*12)
2075 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
2076 t_next = set_date(yearp, climmonth, climday, climhour, climminute, climsecond)
2077 call time_interp(t_next, month, clim_type%FMS_INTP_TYPE_%tweight2, taum, taup, err_msg=err_msg ) !tweight1 is
2078 !! the time weight between the climatology years.
2079 if (taum==2 .and. taup==2) &
2080 clim_type%FMS_INTP_TYPE_%tweight2 = 1._lkind ! protect against post-perth time_interp behavior
2081 if(trim(err_msg) /= '') then
2082 call mpp_error(fatal,'interpolator_3D 5: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2083 endif
2084
2085
2086 if (indexm == clim_type%indexm(i) .and. &
2087 indexp == clim_type%indexp(i) .and. &
2088 climatology == clim_type%climatology(i)) then
2089 else
2090 clim_type%indexm(i) = indexm
2091 clim_type%indexp(i) = indexp
2092 clim_type%climatology(i) = climatology
2093 call read_data(clim_type,clim_type%field_name(i), &
2094 clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), &
2095 clim_type%indexm(i)+(clim_type%climatology(i)-1)*12,i,time)
2096! Read the data for the next month in the previous climatology.
2097 call read_data(clim_type,clim_type%field_name(i), &
2098 clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), &
2099 clim_type%indexp(i)+(clim_type%climatology(i)-1)*12,i,time)
2100 call read_data(clim_type,clim_type%field_name(i), &
2101 clim_type%FMS_INTP_TYPE_%pmon_nyear(:,:,:,i), &
2102 clim_type%indexm(i)+clim_type%climatology(i)*12,i,time)
2103 call read_data(clim_type,clim_type%field_name(i), &
2104 clim_type%FMS_INTP_TYPE_%nmon_nyear(:,:,:,i), &
2105 clim_type%indexp(i)+clim_type%climatology(i)*12,i,time)
2106 endif
2107
2108
2109
2110
2111 else ! We are within a climatology data set
2112
2113 if (taum /= clim_type%time_init(i,1) .or. &
2114 taup /= clim_type%time_init(i,2) ) then
2115
2116 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), taum,i,time)
2117! Read the data for the next month in the previous climatology.
2118 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), taup,i,time)
2119!RSHbug clim_type%pmon_nyear = 0.0
2120!RSHbug clim_type%nmon_nyear = 0.0
2121
2122! clim_type%pmon_nyear(:,:,:,i) = 0.0
2123! clim_type%nmon_nyear(:,:,:,i) = 0.0
2124
2125! set to zero so when next return to bilinear section will be sure to
2126! have proper data (relevant when running fixed_year case for more than
2127! one year in a single job)
2128 clim_type%indexm(i) = 0
2129 clim_type%indexp(i) = 0
2130 clim_type%climatology(i) = 0
2131
2132
2133 clim_type%time_init(i,1) = taum
2134 clim_type%time_init(i,2) = taup
2135 endif
2136! clim_type%tweight3 = 0.0 ! This makes [pn]mon_nyear irrelevant. Set them to 0 to test.
2137 clim_type%FMS_INTP_TYPE_%tweight1 = 0.0_lkind ; clim_type%FMS_INTP_TYPE_%tweight2 = 0.0_lkind
2138 clim_type%FMS_INTP_TYPE_%tweight3 = clim_type%FMS_INTP_TYPE_%tweight
2139 endif
2140
2141 endif ! (BILINEAR)
2142
2143
2144 if(clim_type%TIME_FLAG .eq. linear .and. &
2145 (.not. read_all_on_init) ) then
2146! We need 2 time levels. Check we have the correct data.
2147 clim_type%itaum=0
2148 clim_type%itaup=0
2149 do n=1,size(clim_type%time_init,2)
2150 if (clim_type%time_init(i,n) .eq. taum ) clim_type%itaum = n
2151 if (clim_type%time_init(i,n) .eq. taup ) clim_type%itaup = n
2152 enddo
2153
2154 if (clim_type%itaum.eq.0 .and. clim_type%itaup.eq.0) then
2155!Neither time is set so we need to read 2 time slices.
2156!Set up
2157! field(:,:,:,1) as the previous time slice.
2158! field(:,:,:,2) as the next time slice.
2159 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,1,i), taum,i,time)
2160 clim_type%time_init(i,1) = taum
2161 clim_type%itaum = 1
2162 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,2,i), taup,i,time)
2163 clim_type%time_init(i,2) = taup
2164 clim_type%itaup = 2
2165 endif ! clim_type%itaum.eq.clim_type%itaup.eq.0
2166 if (clim_type%itaum.eq.0 .and. clim_type%itaup.ne.0) then
2167! Can't think of a situation where we would have the next time level but not the previous.
2168 call mpp_error(fatal,'interpolator_3D : No data from the previous climatology time &
2169 & but we have the next time. How did this happen?')
2170 endif
2171 if (clim_type%itaum.ne.0 .and. clim_type%itaup.eq.0) then
2172!We have the previous time step but not the next time step data
2173 clim_type%itaup = 1
2174 if (clim_type%itaum .eq. 1 ) clim_type%itaup = 2
2175 call read_data(clim_type,clim_type%field_name(i), &
2176 clim_type%FMS_INTP_TYPE_%data5d(:,:,:,clim_type%itaup,i), taup,i, time)
2177 clim_type%time_init(i,clim_type%itaup)=taup
2178 endif
2179
2180
2181 endif! TIME_FLAG
2182
2183 endif !( .not. clim_type%separate_time_vary_calc)
2184select case(clim_type%TIME_FLAG)
2185 case (linear)
2186 hinterp_data = (1._lkind-clim_type%FMS_INTP_TYPE_%tweight) * &
2187 clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,clim_type%itaum,i) + &
2188 clim_type%FMS_INTP_TYPE_%tweight * &
2189 clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,clim_type%itaup,i)
2190! case (SEASONAL)
2191! Do sine fit to data at this point
2192 case (bilinear)
2193 hinterp_data = &
2194 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight1)*&
2195 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight3)*clim_type%FMS_INTP_TYPE_%pmon_pyear(istart:iend,jstart:jend,:,i) + &
2196 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight2) * clim_type%FMS_INTP_TYPE_%tweight3 &
2197 * clim_type%FMS_INTP_TYPE_%nmon_pyear(istart:iend,jstart:jend,:,i) + &
2198 clim_type%FMS_INTP_TYPE_%tweight1 * (1._lkind-clim_type%FMS_INTP_TYPE_%tweight3) * &
2199 clim_type%FMS_INTP_TYPE_%pmon_nyear(istart:iend,jstart:jend,:,i) + &
2200 clim_type%FMS_INTP_TYPE_%tweight2 * clim_type%FMS_INTP_TYPE_%tweight3 * &
2201 clim_type%FMS_INTP_TYPE_%nmon_nyear(istart:iend,jstart:jend,:,i)
2202
2203
2204
2205end select
2206
2207select case(clim_type%level_type)
2208 case(pressure)
2209 p_fact = 1.0_lkind
2210 case(sigma)
2211 p_fact = maxval(phalf,3)! max pressure in the column !(:,:,size(phalf,3))
2212end select
2213
2214
2215
2216col_data(:,:)=0.0_lkind
2217select case(clim_type%mr(i))
2218 case(no_conv)
2219 do k = 1,size(hinterp_data,3)
2220 col_data(:,:) = col_data(:,:) + hinterp_data(:,:,k)* &
2221 (clim_type%FMS_INTP_TYPE_%halflevs(k+1)-clim_type%FMS_INTP_TYPE_%halflevs(k))/real(grav,fms_intp_kind_)
2222 enddo
2223
2224 case(kg_m2)
2225 do k = 1,size(hinterp_data,3)
2226 col_data(:,:) = col_data(:,:) + hinterp_data(:,:,k)
2227 hinterp_data(:,:,k) = hinterp_data(:,:,k)/ &
2228 ((clim_type%FMS_INTP_TYPE_%halflevs(k+1)-clim_type%FMS_INTP_TYPE_%halflevs(k))*p_fact)
2229 enddo
2230end select
2231
2232found = .false.
2233do j = 1,size(climo_diag_name(:))
2234 if ( trim(adjustl(lowercase(climo_diag_name(j)))) .eq. trim(adjustl(lowercase(clim_type%field_name(i))))) then
2235 found = .true.
2236 exit
2237 endif
2238enddo
2239
2240if (found) then
2241 if (hinterp_id(j) > 0 ) then
2242 result = send_data(hinterp_id(j),col_data,time,is_in=istart,js_in=jstart)
2243 endif
2244endif
2245
2246
2247!++lwh
2248do j = 1, size(phalf,2)
2249 do ilon=1,size(phalf,1)
2250 pclim = p_fact(ilon,j)*clim_type%FMS_INTP_TYPE_%halflevs
2251 if ( maxval(phalf(ilon,j,:)) > maxval(pclim) ) then
2252 if (verbose > 3) then
2253 call mpp_error(note,"Interpolator: model surface pressure&
2254 & is greater than climatology surface pressure for "&
2255 // trim(clim_type%file_name))
2256 endif
2257 select case(clim_type%out_of_bounds(i))
2258 case(constant)
2259 pclim( maxloc(pclim) ) = maxval( phalf(ilon,j,:) )
2260! case(ZERO)
2261! pclim( maxloc(pclim)) = 0
2262 end select
2263 endif
2264 if ( minval(phalf(ilon,j,:)) < minval(pclim) ) then
2265 if (verbose > 3) then
2266 call mpp_error(note,"Interpolator: model top pressure&
2267 & is less than climatology top pressure for "&
2268 // trim(clim_type%file_name))
2269 endif
2270 select case(clim_type%out_of_bounds(i))
2271 case(constant)
2272 pclim( minloc(pclim) ) = minval( phalf(ilon,j,:) )
2273! case(ZERO)
2274! pclim( maxloc(pclim)) = 0
2275 end select
2276 endif
2277 select case(clim_type%vert_interp(i))
2278 case(interp_weighted_p)
2279 call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:))
2280 case(interp_linear_p)
2281 call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:))
2282! case(INTERP_LOG)
2283 end select
2284 enddo
2285enddo
2286!--lwh
2287
2288select case(clim_type%mr(i))
2289 case(kg_m2)
2290 do k = 1,size(interp_data,3)
2291 interp_data(:,:,k) = interp_data(:,:,k)*(phalf(:,:,k+1)-phalf(:,:,k))
2292 enddo
2293end select
2294 endif !field_name
2295enddo !End of i loop
2296if( .not. found_field) then !field name is not in interpolator file.ERROR.
2297 call mpp_error(fatal,"Interpolator: the field name is not contained in this &
2298 &intepolate_type: "//trim(field_name))
2299endif
2300end subroutine interpolator_3d_
2301!
2302!#######################################################################
2303!
2304!++lwh
2305!---------------------------------------------------------------------
2306!> @brief interpolator_2D receives a field name as input and
2307!! interpolates the field to model a 2D grid and time axis.
2308!!
2309!! @param [inout] <clim_type> The interpolate type previously defined by a call to interpolator_init
2310!! @param [in] <field_name> The name of a field that you wish to interpolate
2311!! @param [in] <Time> The model time that you wish to interpolate to
2312!! @param [in] <is> OPTIONAL: Index for the physics window
2313!! @param [in] <js> OPTIONAL: Index for the physics window
2314!! @param [out] <interp_data> The model fields with the interpolated climatology data
2315!! @param [out] <clim_units> OPTIONAL: The units of field_name
2316!!
2317!! @throw FATAL "interpolator_2D : You must call interpolator_init
2318!! before calling interpolator"
2319!! @throw FATAL "interpolator_2D 1: file="
2320!! @throw FATAL "interpolator_2D 2: file="
2321!! @throw FATAL "interpolator_2D 3: file="
2322!! @throw FATAL "interpolator_2D 4: file="
2323!! @throw FATAL "interpolator_2D 5: file="
2324!! @throw FATAL "interpolator_2D : No data from the previous climatology
2325!! time but we have the next time. How did this happen?"
2326!! @throw FATAL "Interpolator: the field name is not contained in this
2327!! intepolate_type: "
2328subroutine interpolator_2d_(clim_type, Time, interp_data, field_name, is, js, clim_units)
2329!
2330! Return 2-D field interpolated to model grid and time
2331!
2332!
2333! INTENT INOUT
2334! clim_type : The interpolate type previously defined by a call to interpolator_init
2335!
2336! INTENT IN
2337! field_name : The name of the field that you wish to interpolate.
2338! Time : The model time that you wish to interpolate to.
2339! is, js : The indices of the physics window.
2340!
2341! INTENT OUT
2342! interp_data : The model field with the interpolated climatology data.
2343! clim_units : The units of field_name
2344!
2345type(interpolate_type), intent(inout) :: clim_type
2346character(len=*) , intent(in) :: field_name
2347type(time_type) , intent(in) :: Time
2348real(FMS_INTP_KIND_), dimension(:,:), intent(out) :: interp_data
2349integer , intent(in) , optional :: is,js
2350character(len=*) , intent(out), optional :: clim_units
2351integer :: taum, taup
2352real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:)))
2353integer :: istart,iend,jstart,jend
2354logical :: result, found
2355logical :: found_field=.false.
2356integer :: modyear, modmonth, modday, modhour, modminute, modsecond
2357integer :: climyear, climmonth, climday, climhour, climminute, climsecond
2358integer :: year1, month1, day, hour, minute, second
2359integer :: climatology, m
2360type(time_type) :: t_prev, t_next
2361type(time_type), dimension(2) :: month
2362integer :: indexm, indexp, yearm, yearp
2363integer :: j, i, n
2364character(len=256) :: err_msg
2365
2366integer, parameter :: lkind=fms_intp_kind_
2367
2368if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) &
2369 call mpp_error(fatal, "interpolator_2D : You must call interpolator_init before calling interpolator")
2370
2371istart = 1
2372if (present(is)) istart = is
2373iend = istart - 1 + size(interp_data,1)
2374
2375jstart = 1
2376if (present(js)) jstart = js
2377jend = jstart - 1 + size(interp_data,2)
2378
2379do i= 1,size(clim_type%field_name(:))
2380!++lwh
2381 if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then
2382!--lwh
2383
2384 found_field=.true.
2385
2386 if(present(clim_units)) then
2387 call get_variable_units(clim_type%fileobj, clim_type%field_name(i), clim_units)
2388 clim_units = chomp(clim_units)
2389 endif
2390
2391!----------------------------------------------------------------------
2392! skip the time interpolation portion of this routine if subroutine
2393! obtain_interpolator_time_slices has already been called on this
2394! stewp for this interpolate_type variable.
2395!----------------------------------------------------------------------
2396
2397if ( .not. clim_type%separate_time_vary_calc) then
2398! print *, 'TIME INTERPOLATION NOT SEPARATED 2d--', &
2399! trim(clim_type%file_name), mpp_pe()
2400 if (clim_type%climatological_year) then
2401!++lwh
2402 if (size(clim_type%time_slice) > 1) then
2403 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, &
2404 taum, taup, modtime=year, err_msg=err_msg )
2405 if(trim(err_msg) /= '') then
2406 call mpp_error(fatal,'interpolator_2D 1: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2407 endif
2408 else
2409 taum = 1
2410 taup = 1
2411 clim_type%FMS_INTP_TYPE_%tweight = 0._lkind
2412 end if
2413!--lwh
2414 else
2415 call time_interp(time, clim_type%time_slice, clim_type%FMS_INTP_TYPE_%tweight, taum, taup, err_msg=err_msg )
2416 if(trim(err_msg) /= '') then
2417 call mpp_error(fatal,'interpolator_2D 2: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2418 endif
2419 endif
2420
2421! If the climatology file has seasonal, a split time-line or has all the data
2422! read in then enter this loop.
2423!
2424 if(clim_type%TIME_FLAG .ne. linear .or. read_all_on_init) then
2425 clim_type%itaum=taum
2426 clim_type%itaup=taup
2427 endif
2428
2429! if(clim_type%TIME_FLAG .eq. BILINEAR ) then
2430! ! Check if delta-time is greater than delta of first two climatology time-slices.
2431! if ( (Time - clim_type%time_slice(taum) ) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) .or. &
2432! (clim_type%time_slice(taup) - Time) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) ) then
2433! ! The difference between the model time and the last climatology time-slice previous to the model time.
2434! ! We need 2 time levels. Check we have the correct data.
2435! itaum=0
2436! itaup=0
2437! ! Assume this is monthly data. So we need to get the data applicable to the model date but substitute
2438! ! the climatology year into the appropriate place.
2439!
2440! call get_date(Time, modyear, modmonth, modday, modhour, modminute, modsecond)
2441! call get_date(clim_type%time_slice(taum), climyear, climmonth, climday, climhour, climminute, climsecond)
2442! clim_datem = set_date(climyear, modmonth, modday, modhour, modminute, modsecond)
2443! call time_interp(clim_datem, clim_type%time_slice, tweight1, taum1, taup1 )
2444!
2445!
2446! call get_date(clim_type%time_slice(taup), climyear, climmonth, climday, climhour, climminute, climsecond)
2447! clim_datep = set_date(climyear, modmonth, modday, modhour, modminute, modsecond)
2448!
2449!
2450! endif
2451!
2452! endif
2453 if(clim_type%TIME_FLAG .eq. bilinear ) then
2454 ! Check if delta-time is greater than delta of first two climatology time-slices.
2455 if ( (time - clim_type%time_slice(taum) ) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) .or. &
2456 (clim_type%time_slice(taup) - time) > ( clim_type%time_slice(2)- clim_type%time_slice(1) ) ) then
2457 ! The difference between the model time and the last climatology time-slice previous to the model time.
2458 ! We need 2 time levels.
2459 clim_type%itaum=0
2460 clim_type%itaup=0
2461 ! Assume this is monthly data. So we need to get the data applicable to the model date but substitute
2462 ! the climatology year into the appropriate place.
2463
2464
2465 ! We need to get the previous months data for the climatology year before
2466 ! and after the model year.
2467 call get_date(time, modyear, modmonth, modday, modhour, modminute, modsecond)
2468 call get_date(clim_type%time_slice(taum), climyear, climmonth, climday, climhour, climminute, climsecond)
2469
2470 climatology = 1
2471 do m = 1, size(clim_type%clim_times(:,:),2)
2472 !Assume here that a climatology is for 1 year and consists of 12 months starting in January.
2473 call get_date(clim_type%clim_times(1,m), year1, month1, day, hour, minute, second)
2474 if (year1 == climyear) climatology = m
2475 enddo
2476 do m = 1,12
2477 !Find which month we are trying to look at and set clim_date[mp] to the dates spanning that.
2478 call get_date(clim_type%clim_times(m,climatology), year1, month1, day, hour, minute, second)
2479 if ( month1 == modmonth ) then
2480!RSHBUGFX if ( modday <= day ) then
2481 if ( modday < day ) then
2482 indexm = m-1 ; indexp = m
2483 else
2484 indexm = m ; indexp = m+1
2485 endif
2486 endif
2487
2488 enddo
2489 if ( indexm == 0 ) then
2490 indexm = 12
2491 yearm = modyear - 1
2492 else
2493 yearm = modyear
2494 endif
2495 call get_date(clim_type%time_slice(indexm+(climatology-1)*12), &
2496 climyear, climmonth, climday, climhour, climminute, climsecond)
2497 month(1) = set_date(yearm, indexm, climday, climhour, climminute, climsecond)
2498 if ( indexp == 13 ) then
2499 indexp = 1
2500 yearp = modyear + 1
2501 else
2502 yearp = modyear
2503 endif
2504 call get_date(clim_type%time_slice(indexp+(climatology-1)*12), &
2505 climyear, climmonth, climday, climhour, climminute, climsecond)
2506 month(2) = set_date(yearp, indexp, climday, climhour, climminute, climsecond)
2507
2508 call time_interp(time, month, clim_type%FMS_INTP_TYPE_%tweight3, taum, taup, err_msg=err_msg ) ! tweight3 is
2509 !! the time weight between the months.
2510 if (taum==2 .and. taup==2) &
2511 clim_type%FMS_INTP_TYPE_%tweight3 = 1._lkind ! protect against post-perth time_interp behavior
2512 if(trim(err_msg) /= '') then
2513 call mpp_error(fatal,'interpolator_2D 3: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2514 endif
2515
2516 month(1) = clim_type%time_slice(indexm+(climatology-1)*12)
2517 month(2) = clim_type%time_slice(indexm+climatology*12)
2518 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
2519 t_prev = set_date(yearm, climmonth, climday, climhour, climminute, climsecond)
2520 call time_interp(t_prev, month, clim_type%FMS_INTP_TYPE_%tweight1, taum, taup, err_msg=err_msg ) !tweight1 is
2521 !! the time weight between the climatology years.
2522 !> protect against post-perth time_interp behavior below
2523 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight1 = 1._lkind
2524 if(trim(err_msg) /= '') then
2525 call mpp_error(fatal,'interpolator_2D 4: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2526 endif
2527
2528 month(1) = clim_type%time_slice(indexp+(climatology-1)*12)
2529 month(2) = clim_type%time_slice(indexp+climatology*12)
2530 call get_date(month(1), climyear, climmonth, climday, climhour, climminute, climsecond)
2531 t_next = set_date(yearp, climmonth, climday, climhour, climminute, climsecond)
2532 call time_interp(t_next, month, clim_type%FMS_INTP_TYPE_%tweight2, taum, taup, err_msg=err_msg ) !tweight1 is
2533 !! the time weight between the climatology years.
2534 !> protect against post-perth time_interp behavior below
2535 if (taum==2 .and. taup==2) clim_type%FMS_INTP_TYPE_%tweight2 = 1._lkind
2536 if(trim(err_msg) /= '') then
2537 call mpp_error(fatal,'interpolator_2D 5: '//trim(err_msg)//' file='//trim(clim_type%file_name), fatal)
2538 endif
2539
2540
2541 if (indexm == clim_type%indexm(i) .and. &
2542 indexp == clim_type%indexp(i) .and. &
2543 climatology == clim_type%climatology(i)) then
2544 else
2545 clim_type%indexm(i) = indexm
2546 clim_type%indexp(i) = indexp
2547 clim_type%climatology(i) = climatology
2548 call read_data(clim_type,clim_type%field_name(i), &
2549 clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), &
2550 clim_type%indexm(i)+(clim_type%climatology(i)-1)*12,i,time)
2551! Read the data for the next month in the previous climatology.
2552 call read_data(clim_type,clim_type%field_name(i), &
2553 clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), &
2554 clim_type%indexp(i)+(clim_type%climatology(i)-1)*12,i,time)
2555 call read_data(clim_type,clim_type%field_name(i), &
2556 clim_type%FMS_INTP_TYPE_%pmon_nyear(:,:,:,i), &
2557 clim_type%indexm(i)+clim_type%climatology(i)*12,i,time)
2558 call read_data(clim_type,clim_type%field_name(i), &
2559 clim_type%FMS_INTP_TYPE_%nmon_nyear(:,:,:,i), &
2560 clim_type%indexp(i)+clim_type%climatology(i)*12,i,time)
2561 endif
2562
2563
2564
2565
2566 else ! We are within a climatology data set
2567
2568 if (taum /= clim_type%time_init(i,1) .or. &
2569 taup /= clim_type%time_init(i,2) ) then
2570
2571 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), taum,i,time)
2572! Read the data for the next month in the previous climatology.
2573 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%nmon_pyear(:,:,:,i), taup,i,time)
2574!RSHbug clim_type%pmon_nyear = 0.0
2575!RSHbug clim_type%nmon_nyear = 0.0
2576
2577! clim_type%pmon_nyear(:,:,:,i) = 0.0
2578! clim_type%nmon_nyear(:,:,:,i) = 0.0
2579
2580! set to zero so when next return to bilinear section will be sure to
2581! have proper data (relevant when running fixed_year case for more than
2582! one year in a single job)
2583 clim_type%indexm(i) = 0
2584 clim_type%indexp(i) = 0
2585 clim_type%climatology(i) = 0
2586
2587
2588 clim_type%time_init(i,1) = taum
2589 clim_type%time_init(i,2) = taup
2590 endif
2591! clim_type%tweight3 = 0.0 ! This makes [pn]mon_nyear irrelevant. Set them to 0 to test.
2592 clim_type%FMS_INTP_TYPE_%tweight1 = 0.0_lkind
2593 clim_type%FMS_INTP_TYPE_%tweight2 = 0.0_lkind
2594 clim_type%FMS_INTP_TYPE_%tweight3 = clim_type%FMS_INTP_TYPE_%tweight
2595 endif
2596
2597 endif ! (BILINEAR)
2598
2599 if(clim_type%TIME_FLAG .eq. linear .and. &
2600 (.not. read_all_on_init) ) then
2601! We need 2 time levels. Check we have the correct data.
2602 clim_type%itaum=0
2603 clim_type%itaup=0
2604 do n=1,size(clim_type%time_init,2)
2605 if (clim_type%time_init(i,n) .eq. taum ) clim_type%itaum = n
2606 if (clim_type%time_init(i,n) .eq. taup ) clim_type%itaup = n
2607 enddo
2608
2609 if (clim_type%itaum.eq.0 .and. clim_type%itaup.eq.0) then
2610 !Neither time is set so we need to read 2 time slices.
2611 !Set up
2612 ! field(:,:,:,1) as the previous time slice.
2613 ! field(:,:,:,2) as the next time slice.
2614 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,1,i), taum,i,time)
2615 clim_type%time_init(i,1) = taum
2616 clim_type%itaum = 1
2617 call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,2,i), taup,i,time)
2618 clim_type%time_init(i,2) = taup
2619 clim_type%itaup = 2
2620 endif ! clim_type%itaum.eq.clim_type%itaup.eq.0
2621 if (clim_type%itaum.eq.0 .and. clim_type%itaup.ne.0) then
2622 ! Can't think of a situation where we would have the next time level but not the previous.
2623 call mpp_error(fatal,'interpolator_2D : No data from the previous climatology time but we have&
2624 & the next time. How did this happen?')
2625 endif
2626 if (clim_type%itaum.ne.0 .and. clim_type%itaup.eq.0) then
2627 !We have the previous time step but not the next time step data
2628 clim_type%itaup = 1
2629 if (clim_type%itaum .eq. 1 ) clim_type%itaup = 2
2630 call read_data(clim_type,clim_type%field_name(i), &
2631 clim_type%FMS_INTP_TYPE_%data5d(:,:,:,clim_type%itaup,i), taup,i, time)
2632 clim_type%time_init(i,clim_type%itaup)=taup
2633 endif
2634 endif! TIME_FLAG .eq. LINEAR .and. (.not. read_all_on_init)
2635
2636 endif ! (.not. separate_time_vary_calc)
2637
2638
2639
2640select case(clim_type%TIME_FLAG)
2641 case (linear)
2642 hinterp_data = (1._lkind-clim_type%FMS_INTP_TYPE_%tweight)*&
2643 clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,clim_type%itaum,i) &
2644 + clim_type%FMS_INTP_TYPE_%tweight*clim_type%FMS_INTP_TYPE_%data5d &
2645 (istart:iend,jstart:jend,:,clim_type%itaup,i)
2646! case (SEASONAL)
2647! Do sine fit to data at this point
2648 case (bilinear)
2649 hinterp_data = &
2650 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight1) * &
2651 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight3) * &
2652 clim_type%FMS_INTP_TYPE_%pmon_pyear(istart:iend,jstart:jend,:,i) + &
2653 (1._lkind-clim_type%FMS_INTP_TYPE_%tweight2) * &
2654 clim_type%FMS_INTP_TYPE_%tweight3 * clim_type%FMS_INTP_TYPE_%nmon_pyear(istart:iend,jstart:jend,:,i) + &
2655 clim_type%FMS_INTP_TYPE_%tweight1 * (1._lkind-clim_type%FMS_INTP_TYPE_%tweight3) * &
2656 clim_type%FMS_INTP_TYPE_%pmon_nyear(istart:iend,jstart:jend,:,i) + &
2657 clim_type%FMS_INTP_TYPE_%tweight2 * clim_type%FMS_INTP_TYPE_%tweight3 * &
2658 clim_type%FMS_INTP_TYPE_%nmon_nyear(istart:iend,jstart:jend,:,i)
2659
2660end select
2661
2662found = .false.
2663do j = 1,size(climo_diag_name(:))
2664 if (trim(adjustl(lowercase(climo_diag_name(j)))) .eq. trim(adjustl(lowercase(clim_type%field_name(i))))) then
2665 found = .true.
2666 exit
2667 endif
2668enddo
2669
2670if (found) then
2671 if (hinterp_id(j) > 0 ) then
2672 result = send_data(hinterp_id(j),hinterp_data,time,is_in=istart,js_in=jstart)
2673 endif
2674endif
2675
2676 interp_data(:,:) = hinterp_data(:,:,1)
2677
2678 endif !field_name
2679enddo !End of i loop
2680
2681if( .not. found_field) then !field name is not in interpolator file.ERROR.
2682 call mpp_error(fatal,"Interpolator: the field name is not contained in this &
2683 &intepolate_type: "//trim(field_name))
2684endif
2685end subroutine interpolator_2d_
2686!--lwh
2687!
2688!#######################################################################
2689!
2690!---------------------------------------------------------------------
2691!> @brief interpolator_4D_no_time_axis receives a field name as input and
2692!! interpolates the field to model a 4D grid.
2693!!
2694!! @param [inout] <clim_type> The interpolate type previously defined by a call to interpolator_init
2695!! @param [in] <field_name> The name of a field that you wish to interpolate
2696!! @param [in] <phalf> The half level model pressure field
2697!! @param [in] <is> OPTIONAL: Index for the physics window
2698!! @param [in] <js> OPTIONAL: Index for the physics window
2699!! @param [out] <interp_data> The model fields with the interpolated climatology data
2700!! @param [out] <clim_units> OPTIONAL: The units of field_name
2701!!
2702!! @throw FATAL "interpolator_4D_no_time_axis : You must call
2703!! interpolator_init before calling interpolator"
2704!! @throw FATAL "interpolator_mod: cannot use 4D interface to
2705!! interpolator for this file"
2706!! @throw NOTE "Interpolator: model surface pressure is greater than
2707!! surface pressure of input data for "
2708!! @throw NOTE "Interpolator: model top pressure is less than surface
2709!! pressure of input data for "
2710!! @throw FATAL "Interpolator: the field name is not contained in this
2711!! intepolate_type: "
2712subroutine interpolator_4d_no_time_axis_(clim_type, phalf, interp_data, field_name, is,js, clim_units)
2713
2714! Return 4-D field interpolated to model grid
2715
2716! INTENT INOUT
2717! clim_type : The interpolate type previously defined by a call to interpolator_init
2718
2719! INTENT IN
2720! field_name : The name of a field that you wish to interpolate.
2721! all variables within this interpolate_type variable
2722! will be interpolated on this call. field_name may
2723! be any one of the variables.
2724! phalf : The half level model pressure field.
2725! is, js : The indices of the physics window.
2726
2727! INTENT OUT
2728! interp_data : The model fields
2729! clim_units : The units of field_name
2730
2731type(interpolate_type), intent(inout) :: clim_type
2732character(len=*) , intent(in) :: field_name
2733real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf
2734real(FMS_INTP_KIND_), dimension(:,:,:,:), intent(out) :: interp_data
2735integer , intent(in) , optional :: is,js
2736character(len=*) , intent(out), optional :: clim_units
2737integer :: ilon
2738real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),&
2739 size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:)),size(clim_type%field_name(:)))
2740real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2))
2741real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:)))
2742integer :: istart,iend,jstart,jend
2743logical :: found_field=.false.
2744integer :: i, j, k, n
2745
2746integer, parameter :: lkind=fms_intp_kind_
2747
2748if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) &
2749 call mpp_error(fatal, "interpolator_4D_no_time_axis : You must call interpolator_init before calling interpolator")
2750
2751do n=2,size(clim_type%field_name(:))
2752 if (clim_type%vert_interp(n) /= clim_type%vert_interp(n-1) .or. &
2753 clim_type%out_of_bounds(n) /= clim_type%out_of_bounds(n-1)) then
2754 if (mpp_pe() == mpp_root_pe() ) then
2755 print *, 'processing file ' // trim(clim_type%file_name)
2756 endif
2757 call mpp_error (fatal, 'interpolator_mod: &
2758 &cannot use 4D interface to interpolator for this file')
2759 endif
2760end do
2761
2762istart = 1
2763if (present(is)) istart = is
2764iend = istart - 1 + size(interp_data,1)
2765
2766jstart = 1
2767if (present(js)) jstart = js
2768jend = jstart - 1 + size(interp_data,2)
2769
2770do i= 1,size(clim_type%field_name(:))
2771 if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then
2772 found_field=.true.
2773 exit
2774 endif
2775end do
2776i = 1
2777
2778if(present(clim_units)) then
2779 call get_variable_units(clim_type%fileobj, clim_type%field_name(i), clim_units)
2780 clim_units = chomp(clim_units)
2781endif
2782
2783do n=1, size(clim_type%field_name(:))
2784 hinterp_data(:,:,:,n) = clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,1,n)
2785end do
2786
2787select case(clim_type%level_type)
2788 case(pressure)
2789 p_fact = 1.0_lkind
2790 case(sigma)
2791 p_fact = maxval(phalf,3)! max pressure in the column !(:,:,size(phalf,3))
2792end select
2793
2794 do i= 1, size(clim_type%field_name(:))
2795 select case(clim_type%mr(i))
2796 case(kg_m2)
2797 do k = 1,size(hinterp_data,3)
2798 hinterp_data(:,:,k,i) = hinterp_data(:,:,k,i)/&
2799 ((clim_type%FMS_INTP_TYPE_%halflevs(k+1)-clim_type%FMS_INTP_TYPE_%halflevs(k))*p_fact)
2800 enddo
2801 end select
2802 enddo
2803
2804 i = 1
2805
2806do j = 1, size(phalf,2)
2807 do ilon=1,size(phalf,1)
2808 pclim = p_fact(ilon,j)*clim_type%FMS_INTP_TYPE_%halflevs
2809 if ( maxval(phalf(ilon,j,:)) > maxval(pclim) ) then
2810 if (verbose > 3) then
2811 call mpp_error(note,"Interpolator: model surface pressure&
2812 & is greater than surface pressure of input data for "&
2813 // trim(clim_type%file_name))
2814 endif
2815 select case(clim_type%out_of_bounds(i))
2816 case(constant)
2817 pclim( maxloc(pclim) ) = maxval( phalf(ilon,j,:) )
2818 end select
2819 endif
2820 if ( minval(phalf(ilon,j,:)) < minval(pclim) ) then
2821 if (verbose > 3) then
2822 call mpp_error(note,"Interpolator: model top pressure&
2823 & is less than top pressure of input data for "&
2824 // trim(clim_type%file_name))
2825 endif
2826 select case(clim_type%out_of_bounds(i))
2827 case(constant)
2828 pclim( minloc(pclim) ) = minval( phalf(ilon,j,:) )
2829 end select
2830 endif
2831 select case(clim_type%vert_interp(i))
2832 case(interp_weighted_p)
2833 call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,:),interp_data(ilon,j,:,:))
2834 case(interp_linear_p)
2835 do n=1, size(clim_type%field_name(:))
2836 call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,n),interp_data(ilon,j,:,n))
2837 end do
2838 end select
2839 enddo
2840enddo
2841
2842 do i= 1, size(clim_type%field_name(:))
2843
2844select case(clim_type%mr(i))
2845 case(kg_m2)
2846 do k = 1,size(interp_data,3)
2847 interp_data(:,:,k,i) = interp_data(:,:,k,i)*(phalf(:,:,k+1)-phalf(:,:,k))
2848 enddo
2849end select
2850
2851 end do
2852
2853if( .not. found_field) then !field name is not in interpolator file.ERROR.
2854 call mpp_error(fatal,"Interpolator: the field name is not contained in this &
2855 &intepolate_type: "//trim(field_name))
2856endif
2857end subroutine interpolator_4d_no_time_axis_
2858
2859!#######################################################################
2860!
2861!---------------------------------------------------------------------
2862!> @brief interpolator_3D_no_time_axis receives a field name as input and
2863!! interpolates the field to model a 3D grid.
2864!!
2865!! @param [inout] <clim_type> The interpolate type previously defined by a call to interpolator_init
2866!! @param [in] <field_name> The name of a field that you wish to interpolate
2867!! @param [in] <phalf> The half level model pressure field
2868!! @param [in] <is> OPTIONAL: Index for the physics window
2869!! @param [in] <js> OPTIONAL: Index for the physics window
2870!! @param [out] <interp_data> The model fields with the interpolated climatology data
2871!! @param [out] <clim_units> OPTIONAL: The units of field_name
2872!!
2873!! @throw FATAL "interpolator_3D_no_time_axis : You must call
2874!! interpolator_init before calling interpolator"
2875!! @throw FATAL "interpolator_mod: cannot use 4D interface to
2876!! interpolator for this file"
2877!! @throw NOTE "Interpolator: model surface pressure is greater than
2878!! climatology surface pressure for "
2879!! @throw NOTE "Interpolator: model top pressure is less than
2880!! climatology top pressure for "
2881!! @throw FATAL "Interpolator: the field name is not contained in this
2882!! intepolate_type: "
2883subroutine interpolator_3d_no_time_axis_(clim_type, phalf, interp_data, field_name, is,js, clim_units)
2884
2885! Return 3-D field interpolated to model grid
2886
2887! INTENT INOUT
2888! clim_type : The interpolate type previously defined by a call to interpolator_init
2889
2890! INTENT IN
2891! field_name : The name of the field that you wish to interpolate.
2892! phalf : The half level model pressure field.
2893! is, js : The indices of the physics window.
2894
2895! INTENT OUT
2896! interp_data : The model field with the interpolated climatology data.
2897! clim_units : The units of field_name
2898
2899type(interpolate_type), intent(inout) :: clim_type
2900character(len=*) , intent(in) :: field_name
2901real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf
2902real(FMS_INTP_KIND_), dimension(:,:,:), intent(out) :: interp_data
2903integer , intent(in) , optional :: is,js
2904character(len=*) , intent(out), optional :: clim_units
2905integer :: ilon !< No description
2906real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),&
2907 size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:))) !< No description
2908real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2)) !< No description
2909real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:))) !< No description
2910integer :: istart,iend,jstart,jend !< No description
2911logical :: found_field=.false. !< No description
2912integer :: i, j, k !< No description
2913
2914integer, parameter :: lkind=fms_intp_kind_
2915
2916if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) &
2917 call mpp_error(fatal, "interpolator_3D_no_time_axis : You must call interpolator_init before calling interpolator")
2918
2919istart = 1
2920if (present(is)) istart = is
2921iend = istart - 1 + size(interp_data,1)
2922
2923jstart = 1
2924if (present(js)) jstart = js
2925jend = jstart - 1 + size(interp_data,2)
2926
2927do i= 1,size(clim_type%field_name(:))
2928 if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then
2929 found_field=.true.
2930 if(present(clim_units)) then
2931 call get_variable_units(clim_type%fileobj, clim_type%field_name(i), clim_units)
2932 clim_units = chomp(clim_units)
2933 endif
2934
2935 hinterp_data = clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,1,i)
2936
2937select case(clim_type%level_type)
2938 case(pressure)
2939 p_fact = 1.0_lkind
2940 case(sigma)
2941 p_fact = maxval(phalf,3)! max pressure in the column !(:,:,size(phalf,3))
2942end select
2943
2944select case(clim_type%mr(i))
2945 case(kg_m2)
2946 do k = 1,size(hinterp_data,3)
2947 hinterp_data(:,:,k) = hinterp_data(:,:,k)/&
2948 ((clim_type%FMS_INTP_TYPE_%halflevs(k+1)-clim_type%FMS_INTP_TYPE_%halflevs(k))*p_fact)
2949 enddo
2950end select
2951
2952do j = 1, size(phalf,2)
2953 do ilon=1,size(phalf,1)
2954 pclim = p_fact(ilon,j)*clim_type%FMS_INTP_TYPE_%halflevs
2955 if ( maxval(phalf(ilon,j,:)) > maxval(pclim) ) then
2956 if (verbose > 3) then
2957 call mpp_error(note,"Interpolator: model surface pressure&
2958 & is greater than climatology surface pressure for "&
2959 // trim(clim_type%file_name))
2960 endif
2961 select case(clim_type%out_of_bounds(i))
2962 case(constant)
2963 pclim( maxloc(pclim) ) = maxval( phalf(ilon,j,:) )
2964 end select
2965 endif
2966 if ( minval(phalf(ilon,j,:)) < minval(pclim) ) then
2967 if (verbose > 3) then
2968 call mpp_error(note,"Interpolator: model top pressure&
2969 & is less than climatology top pressure for "&
2970 // trim(clim_type%file_name))
2971 endif
2972 select case(clim_type%out_of_bounds(i))
2973 case(constant)
2974 pclim( minloc(pclim) ) = minval( phalf(ilon,j,:) )
2975 end select
2976 endif
2977 select case(clim_type%vert_interp(i))
2978 case(interp_weighted_p)
2979 call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:))
2980 case(interp_linear_p)
2981 call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:))
2982 end select
2983 enddo
2984enddo
2985
2986select case(clim_type%mr(i))
2987 case(kg_m2)
2988 do k = 1,size(interp_data,3)
2989 interp_data(:,:,k) = interp_data(:,:,k)*(phalf(:,:,k+1)-phalf(:,:,k))
2990 enddo
2991end select
2992
2993 endif !field_name
2994enddo !End of i loop
2995if( .not. found_field) then !field name is not in interpolator file.ERROR.
2996 call mpp_error(fatal,"Interpolator: the field name is not contained in this &
2997 &intepolate_type: "//trim(field_name))
2998endif
2999end subroutine interpolator_3d_no_time_axis_
3000
3001!#######################################################################
3002!
3003!---------------------------------------------------------------------
3004!> @brief interpolator_2D_no_time_axis receives a field name as input and
3005!! interpolates the field to model a 2D grid.
3006!!
3007!! @param [inout] <clim_type> The interpolate type previously defined by a call to interpolator_init
3008!! @param [in] <field_name> The name of a field that you wish to interpolate
3009!! @param [in] <is> OPTIONAL: Index for the physics window
3010!! @param [in] <js> OPTIONAL: Index for the physics window
3011!! @param [out] <interp_data> The model fields with the interpolated climatology data
3012!! @param [out] <clim_units> OPTIONAL: The units of field_name
3013!!
3014!! @throw FATAL "interpolator_2D_no_time_axis : You must call
3015!! interpolator_init before calling interpolator"
3016!! @throw FATAL "Interpolator: the field name is not contained in this
3017!! intepolate_type: "
3018subroutine interpolator_2d_no_time_axis_(clim_type, interp_data, field_name, is, js, clim_units)
3019
3020! Return 2-D field interpolated to model grid
3021
3022! INTENT INOUT
3023! clim_type : The interpolate type previously defined by a call to interpolator_init
3024
3025! INTENT IN
3026! field_name : The name of the field that you wish to interpolate.
3027! is, js : The indices of the physics window.
3028
3029! INTENT OUT
3030! interp_data : The model field with the interpolated climatology data.
3031! clim_units : The units of field_name
3032
3033type(interpolate_type), intent(inout) :: clim_type
3034character(len=*) , intent(in) :: field_name
3035real(FMS_INTP_KIND_), dimension(:,:), intent(out) :: interp_data
3036integer , intent(in) , optional :: is,js
3037character(len=*) , intent(out), optional :: clim_units
3038real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),&
3039 size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:)))
3040integer :: istart,iend,jstart,jend
3041logical :: found_field=.false.
3042integer :: i
3043
3044if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) &
3045 call mpp_error(fatal, "interpolator_2D_no_time_axis : You must call interpolator_init before calling interpolator")
3046
3047istart = 1
3048if (present(is)) istart = is
3049iend = istart - 1 + size(interp_data,1)
3050
3051jstart = 1
3052if (present(js)) jstart = js
3053jend = jstart - 1 + size(interp_data,2)
3054
3055do i= 1,size(clim_type%field_name(:))
3056 if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then
3057
3058 found_field=.true.
3059
3060 if(present(clim_units)) then
3061 call get_variable_units(clim_type%fileobj, clim_type%field_name(i), clim_units)
3062 clim_units = chomp(clim_units)
3063 endif
3064
3065 hinterp_data = clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,1,i)
3066
3067 interp_data(:,:) = hinterp_data(:,:,1)
3068
3069 endif !field_name
3070enddo !End of i loop
3071
3072if( .not. found_field) then !field name is not in interpolator file.ERROR.
3073 call mpp_error(fatal,"Interpolator: the field name is not contained in this &
3074 &intepolate_type: "//trim(field_name))
3075endif
3076
3077end subroutine interpolator_2d_no_time_axis_
3078
3079!#######################################################################
3080!
3081!---------------------------------------------------------------------
3082!> @brief read_data receives various climate data as inputs and
3083!! returns a horizontally interpolated climatology field.
3084!!
3085!! @param [in] <clim_type> The interpolate type which contains the data
3086!! @param [in] <src_field> The field type
3087!! @param [in] <nt> The index of the time slice of the climatology that you wish to read
3088!! @param [in] <i> OPTIONAL: The index of the field name that you are trying to read
3089!! @param [in] <Time> OPTIONAL: The model time. Used for diagnostic purposes only
3090!! @param [out] <hdata> The horizontally interpolated climatology field. This
3091! field will still be on the climatology vertical grid
3092subroutine read_data_(clim_type,field_name, hdata, nt,i, Time)
3093!
3094! INTENT IN
3095! clim_type : The interpolate type which contains the data
3096! src_field : The field type
3097! nt : The index of the time slice of the climatology that you wish to read.
3098! i : The index of the field name that you are trying to read. (optional)
3099! Time : The model time. Used for diagnostic purposes only. (optional)
3100!
3101! INTENT OUT
3102!
3103! hdata : The horizontally interpolated climatology field. This
3104! field will still be on the climatology vertical grid.
3105!
3106type(interpolate_type) , intent(in) :: clim_type
3107character(len=*) , intent(in) :: field_name
3108integer , intent(in) :: nt
3109real(FMS_INTP_KIND_) , intent(out) :: hdata(:,:,:)
3110integer , intent(in) :: i
3111type(time_type), optional, intent(in) :: Time
3112
3113integer :: k, km
3114! sjs
3115real(FMS_INTP_KIND_), allocatable :: climdata(:,:,:), climdata2(:,:,:)
3116
3117 allocate(climdata(size(clim_type%FMS_INTP_TYPE_%lon(:)),size(clim_type%FMS_INTP_TYPE_%lat(:)), &
3118 size(clim_type%FMS_INTP_TYPE_%levs(:))))
3119 if(clim_type%has_level(i)) then ! has vertical level
3120 call fms2_io_read_data(clim_type%fileobj,field_name, climdata,nt)
3121 else ! no vertical level
3122 call fms2_io_read_data(clim_type%fileobj,field_name, climdata(:,:,1),nt)
3123 endif
3124
3125! if vertical index increases upward, flip the data so that lowest
3126! pressure level data is at index 1, rather than the highest pressure
3127! level data. the indices themselves were previously flipped.
3128 if (clim_type%vertical_indices == increasing_upward) then
3129 allocate(climdata2(size(clim_type%FMS_INTP_TYPE_%lon(:)), &
3130 size(clim_type%FMS_INTP_TYPE_%lat(:)), &
3131 size(clim_type%FMS_INTP_TYPE_%levs(:))))
3132 km = size(clim_type%FMS_INTP_TYPE_%levs(:))
3133 do k=1, km
3134 climdata2(:,:,k) = climdata(:,:,km+1-k)
3135 end do
3136 climdata = climdata2
3137 deallocate (climdata2)
3138 endif
3139 call horiz_interp(clim_type%interph, climdata, hdata)
3140 if (clim_diag_initialized) &
3141 call diag_read_data(clim_type,climdata,i, time)
3142 deallocate(climdata)
3143
3144
3145 end subroutine read_data_
3146
3147!#######################################################################
3148!
3149!---------------------------------------------------------------------
3150!> @brief read_data_no_time_axis receives various climate data as inputs and
3151!! returns a horizontally interpolated climatology field without the
3152!! time axis.
3153!!
3154!! @param [in] <clim_type> The interpolate type which contains the data
3155!! @param [in] <src_field> The field type
3156!! @param [in] <i> OPTIONAL: The index of the field name that you are trying to read
3157!! @param [out] <hdata> The horizontally interpolated climatology field. This
3158! field will still be on the climatology vertical grid
3159subroutine read_data_no_time_axis_(clim_type,field_name, hdata, i)
3160
3161! INTENT IN
3162! clim_type : The interpolate type which contains the data
3163! src_field : The field type
3164! i : The index of the field name that you are trying to read. (optional)
3165
3166! INTENT OUT
3167
3168! hdata : The horizontally interpolated climatology field. This
3169! field will still be on the climatology vertical grid.
3170
3171type(interpolate_type) , intent(in) :: clim_type
3172character(len=*) , intent(in) :: field_name
3173real(FMS_INTP_KIND_) , intent(out) :: hdata(:,:,:)
3174integer , intent(in) :: i
3175
3176integer :: k, km
3177! sjs
3178real(FMS_INTP_KIND_), allocatable :: climdata(:,:,:), climdata2(:,:,:)
3179
3180 allocate(climdata(size(clim_type%FMS_INTP_TYPE_%lon(:)),&
3181 size(clim_type%FMS_INTP_TYPE_%lat(:)),&
3182 size(clim_type%FMS_INTP_TYPE_%levs(:))))
3183
3184 if(clim_type%has_level(i)) then ! has vertical level
3185 call fms2_io_read_data(clim_type%fileobj,field_name, climdata)
3186 else ! no vertical level
3187 call fms2_io_read_data(clim_type%fileobj,field_name, climdata(:,:,1))
3188 endif
3189! if vertical index increases upward, flip the data so that lowest
3190! pressure level data is at index 1, rather than the highest pressure
3191! level data. the indices themselves were previously flipped.
3192 if (clim_type%vertical_indices == increasing_upward) then
3193 allocate(climdata2(size(clim_type%FMS_INTP_TYPE_%lon(:)), &
3194 size(clim_type%FMS_INTP_TYPE_%lat(:)), &
3195 size(clim_type%FMS_INTP_TYPE_%levs(:))))
3196 km = size(clim_type%FMS_INTP_TYPE_%levs(:))
3197 do k=1, km
3198 climdata2(:,:,k) = climdata(:,:,km+1-k)
3199 end do
3200 climdata = climdata2
3201 deallocate (climdata2)
3202 endif
3203
3204 call horiz_interp(clim_type%interph, climdata, hdata)
3205 deallocate(climdata)
3206
3207end subroutine read_data_no_time_axis_
3208
3209!#######################################################################
3210!
3211!---------------------------------------------------------------------
3212!> @brief diag_read_data receives the data read in by read_data as
3213!! inputs and runs a diagnosis.
3214!!
3215!! @param [in] <clim_type> The interpolate type which contains the data
3216!! @param [in] <model_data> The data read in from file that is being diagnosed
3217!! @param [in] <i> The index of the field name that you are diagnosing
3218!! @param [in] <Time> The model time.
3219subroutine diag_read_data_(clim_type,model_data, i, Time)
3220!
3221! A routine to diagnose the data read in by read_data
3222!
3223! INTENT IN
3224! clim_type : The interpolate type.
3225! model_data : The data read in from file that is being diagnosed.
3226! i : The index of the field name that you are diagnosing.
3227! Time : The model time
3228!
3229type(interpolate_type), intent(in) :: clim_type
3230real(FMS_INTP_KIND_) , intent(in) :: model_data(:,:,:)
3231integer , intent(in) :: i
3232type(time_type) , intent(in) :: Time
3233
3234integer :: j,k
3235real(FMS_INTP_KIND_) :: col_data(size(model_data,1),size(model_data,2))
3236logical :: result, found
3237
3238integer, parameter :: lkind=fms_intp_kind_
3239
3240found = .false.
3241do j = 1,size(climo_diag_name(:))
3242 if (trim(adjustl(lowercase(climo_diag_name(j)))) .eq. trim(adjustl(lowercase(clim_type%field_name(i))))) then
3243 found = .true.
3244 exit
3245 endif
3246enddo
3247
3248if(found) then
3249 if(climo_diag_id(j)>0) then
3250 col_data(:,:)=0.0_lkind
3251 do k=1,size(model_data,3)
3252 col_data(:,:) = col_data(:,:) + &
3253 model_data(:,:,k)* &
3254 (clim_type%FMS_INTP_TYPE_%halflevs(k+1)-clim_type%FMS_INTP_TYPE_%halflevs(k))/real(grav,fms_intp_kind_)
3255 enddo
3256 result = send_data(climo_diag_id(j),col_data(clim_type%is:clim_type%ie,clim_type%js:clim_type%je),time)
3257 endif
3258endif
3259
3260end subroutine diag_read_data_
3261!
3262!#################################################################
3263!
3264!
3265!---------------------------------------------------------------------
3266!> @brief interp_weighted_scalar_2D receives the variables grdin,
3267!! grdout, and datin as inputs and returns datout.
3268!!
3269!! @param [in] <grdin> No description
3270!! @param [in] <grdout> No description
3271!! @param [in] <datin> No description
3272!! @param [out] <datout> No description
3273subroutine interp_weighted_scalar_2d_ (grdin, grdout, datin, datout )
3274real(FMS_INTP_KIND_), intent(in), dimension(:) :: grdin, grdout
3275real(FMS_INTP_KIND_), intent(in), dimension(:,:) :: datin
3276real(FMS_INTP_KIND_), intent(out), dimension(:,:) :: datout
3277
3278integer :: j, k, n
3279integer, parameter :: lkind=fms_intp_kind_
3280
3281if (size(grdin(:)).ne. (size(datin,1)+1)) &
3282 call mpp_error(fatal,'interp_weighted_scalar : input data and pressure do not have the same number of levels')
3283if (size(grdout(:)).ne. (size(datout,1 )+1)) &
3284 call mpp_error(fatal,'interp_weighted_scalar : output data and pressure do not have the same number of levels')
3285
3286 do k = 1, size(datout,1 )
3287 datout(k,:) = 0.0_lkind
3288
3289 do j = 1, size(datin,1 )
3290
3291 if ( grdin(j) <= grdout(k) .and. &
3292 grdin(j+1) >= grdout(k) .and. &
3293 grdin(j+1) <= grdout(k+1) ) then
3294
3295 do n= 1, size(datin,2)
3296 datout(k,n) = datout(k,n) + datin(j,n)*(grdin(j+1)-grdout(k))
3297 end do
3298
3299 else if ( grdin(j) >= grdout(k) .and. &
3300 grdin(j) <= grdout(k+1) .and. &
3301 grdin(j+1) >= grdout(k+1) ) then
3302
3303 do n= 1, size(datin,2)
3304 datout(k,n) = datout(k,n) + datin(j,n)*(grdout(k+1)-grdin(j))
3305 end do
3306
3307 else if ( grdin(j) >= grdout(k) .and. &
3308 grdin(j+1) <= grdout(k+1) ) then
3309
3310 do n= 1, size(datin,2)
3311 datout(k,n) = datout(k,n) + datin(j,n)*(grdin(j+1)-grdin(j))
3312 end do
3313
3314 else if ( grdin(j) <= grdout(k) .and. &
3315 grdin(j+1) >= grdout(k+1) ) then
3316
3317 do n= 1, size(datin,2)
3318 datout(k,n) = datout(k,n) + datin(j,n)*(grdout(k+1)-grdout(k))
3319
3320 end do
3321 endif
3322
3323 enddo
3324
3325 do n= 1, size(datin,2)
3326 datout(k,n) = datout(k,n)/(grdout(k+1)-grdout(k))
3327 end do
3328
3329 enddo
3330
3331end subroutine interp_weighted_scalar_2d_
3332
3333!---------------------------------------------------------------------
3334!> @brief interp_weighted_scalar_1D receives the variables grdin,
3335!! grdout, and datin as inputs and returns datout.
3336!!
3337!! @param [in] <grdin> No description
3338!! @param [in] <grdout> No description
3339!! @param [in] <datin> No description
3340!! @param [out] <datout> No description
3341subroutine interp_weighted_scalar_1d_ (grdin, grdout, datin, datout )
3342real(FMS_INTP_KIND_), intent(in), dimension(:) :: grdin, grdout, datin
3343real(FMS_INTP_KIND_), intent(out), dimension(:) :: datout
3344
3345integer :: j, k
3346integer, parameter :: lkind=fms_intp_kind_
3347
3348if (size(grdin(:)).ne. (size(datin(:))+1)) &
3349 call mpp_error(fatal,'interp_weighted_scalar : input data and pressure do not have the same number of levels')
3350if (size(grdout(:)).ne. (size(datout(:))+1)) &
3351 call mpp_error(fatal,'interp_weighted_scalar : output data and pressure do not have the same number of levels')
3352
3353 do k = 1, size(datout(:))
3354 datout(k) = 0.0_lkind
3355
3356 do j = 1, size(datin(:))
3357
3358 if ( grdin(j) <= grdout(k) .and. &
3359 grdin(j+1) >= grdout(k) .and. &
3360 grdin(j+1) <= grdout(k+1) ) then
3361
3362 datout(k) = datout(k) + datin(j)*(grdin(j+1)-grdout(k))
3363
3364 else if ( grdin(j) >= grdout(k) .and. &
3365 grdin(j) <= grdout(k+1) .and. &
3366 grdin(j+1) >= grdout(k+1) ) then
3367
3368 datout(k) = datout(k) + datin(j)*(grdout(k+1)-grdin(j))
3369
3370 else if ( grdin(j) >= grdout(k) .and. &
3371 grdin(j+1) <= grdout(k+1) ) then
3372
3373 datout(k) = datout(k) + datin(j)*(grdin(j+1)-grdin(j))
3374
3375 else if ( grdin(j) <= grdout(k) .and. &
3376 grdin(j+1) >= grdout(k+1) ) then
3377
3378 datout(k) = datout(k) + datin(j)*(grdout(k+1)-grdout(k))
3379
3380 endif
3381
3382 enddo
3383
3384 datout(k) = datout(k)/(grdout(k+1)-grdout(k))
3385
3386 enddo
3387
3388end subroutine interp_weighted_scalar_1d_
3389!
3390!#################################################################
3391!
3392!---------------------------------------------------------------------
3393!> @brief interp_linear receives the variables grdin,
3394!! grdout, and datin as inputs and returns a linear
3395!! interpolation.
3396!!
3397!! @param [in] <grdin> No description
3398!! @param [in] <grdout> No description
3399!! @param [in] <datin> No description
3400!! @param [out] <datout> No description
3401subroutine interp_linear_ ( grdin, grdout, datin, datout )
3402real(FMS_INTP_KIND_), intent(in), dimension(:) :: grdin, grdout, datin
3403real(FMS_INTP_KIND_), intent(out), dimension(:) :: datout
3404
3405integer :: j, k, n
3406real(FMS_INTP_KIND_) :: wt
3407
3408integer, parameter :: lkind=fms_intp_kind_
3409
3410if (size(grdin(:)).ne. (size(datin(:))+1)) &
3411 call mpp_error(fatal,'interp_linear : input data and pressure do not have the same number of levels')
3412if (size(grdout(:)).ne. (size(datout(:))+1)) &
3413 call mpp_error(fatal,'interp_linear : output data and pressure do not have the same number of levels')
3414
3415
3416 n = size(grdin(:))
3417
3418 do k= 1, size(datout(:))
3419
3420 ! ascending grid values
3421 if (grdin(1) < grdin(n)) then
3422 do j = 2, size(grdin(:))-1
3423 if (grdout(k) <= grdin(j)) exit
3424 enddo
3425 ! descending grid values
3426 else
3427 do j = size(grdin(:)), 3, -1
3428 if (grdout(k) <= grdin(j-1)) exit
3429 enddo
3430 endif
3431
3432 ! linear interpolation
3433 wt = (grdout(k)-grdin(j-1)) / (grdin(j)-grdin(j-1))
3434!print '(a,2i3,4f6.1)', 'k,j=',k,j,grdout(k),grdin(j-1),grdin(j),wt
3435 ! constant value extrapolation
3436 ! wt = min(max(wt,0.),1.)
3437
3438 datout(k) = (1._lkind-wt)*datin(j-1) + wt*datin(j)
3439
3440 enddo
3441
3442end subroutine interp_linear_
3443!
3444!########################################################################
3445
3446!> @}
3447! close documentation grouping
3448!
3449!#######################################################################