!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module forcing_tools 14,9
!BOP
! !MODULE: forcing_tools
! !DESCRIPTION:
! Contains common variables, parameters and routines that are used by
! all of the individual forcing modules.
!
! !REVISION HISTORY:
! SVN:$Id: forcing_tools.F90 808 2006-04-28 17:06:38Z njn01 $
!
! !USES:
use kinds_mod
use blocks
use domain
use constants
use io
use io_types
use grid
use time_management
use exit_mod
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: find_interp_time, &
get_forcing_filename, &
find_forcing_times, &
interpolate_forcing, &
update_forcing_data, &
echo_forcing_options
! !PUBLIC DATA MEMBERS:
real (r8), parameter, public :: &
never = 1.e20_r8, &! value to signify never
always = -1.e20_r8 ! value to signify always
!*** common data shared by forcing routines
!EOP
!BOC
!-----------------------------------------------------------------------
!
! internal module variables
!
!-----------------------------------------------------------------------
real (r8), dimension(12) :: &
thour00_midmonth, &
thour00_endmonth
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: find_interp_time
! !INTERFACE:
subroutine find_interp_time(forcing_interp_inc, forcing_interp_next) 6
! !DESCRIPTION:
! Finds the next time at which temporal interpolation is needed
! in the forcing.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), intent(in) :: &
forcing_interp_inc ! increment (hours) for intepolation time
! !OUTPUT PARAMETERS:
real (r8), intent(out) :: &
forcing_interp_next ! next interpolation time
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: nstart, n
real (r8) :: forcing_time_test
!-----------------------------------------------------------------------
!
! find the first interpolation time that is greater than the current
! time, which then becomes the next interpolation time.
!
!-----------------------------------------------------------------------
nstart = int(thour00/forcing_interp_inc) - 2
do n = nstart, nstart + 5
forcing_time_test = (n + p5)*forcing_interp_inc
if (forcing_time_test > thour00) exit
enddo
forcing_interp_next = forcing_time_test
!-----------------------------------------------------------------------
!EOC
end subroutine find_interp_time
!***********************************************************************
!BOP
! !IROUTINE: get_forcing_filename
! !INTERFACE:
subroutine get_forcing_filename(forcing_filename, forcing_infile, & 7,1
forcing_time, forcing_data_inc)
! !DESCRIPTION:
! Finds the name of an n-hour file given the forcing time needed.
! Files can be labeled by the date of the beginning of the
! forcing interval or by the date of the middle in the forcing
! interval (which corresponds to the value of the variable
! forcing\_time) depending on the value of
! filename\_at\_begin\_of\_interval (T or F).
!
! File names are asumed to be of the form:
! fileroot.year.day.hour (fileroot.yyyy.ddd.hh)
! where fileroot is specified by namelist input, year is
! 4 digits, day is 3 digits, and hour is 2 digits.
! for example, Jan 1 of year 22 with a root of 'ws' would be
! ws.0022.001.00 (labeled by beginning of interval)
! ws.0022.001.12 (labeled by middle of 1-day interval)
!
! Remember, forcing times are relative to 01-01-0000 so results
! may not be what you expect. for example, if the forcing
! increment is 2 days, then even numbered years will have
! mid-interval-referenced files that are on even numbered
! days, and odd numbered years will reference odd numbered
! days (assuming no leap years).
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), intent(in) :: &
forcing_time, &! time at begin of forcing interval
forcing_data_inc ! increment in hours to read data
character (*), intent(in) :: &
forcing_infile ! root filename
! !OUTPUT PARAMETERS:
character (char_len), intent(out) :: &
forcing_filename ! output forcing filename
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
file_year, &! year for file name
file_day, &! day for file name
file_hour, &! hour for file name
n, &! dummy loop index
days_in_prev_year, &! days in previous year
cindx ! index for character strings
real (r8) :: &
file_time, &!
file_time00 !
character (len = 5) :: file_cyear ! character version of the
character (len = 4) :: file_cday ! above integer fields
character (len = 3) :: file_chour
logical (log_kind) :: &
filename_at_begin_of_interval ! flag for position of data file
! in forcing interval (begin,middle)
!-----------------------------------------------------------------------
!
! filename_at_begin_of_interval = .true. for filenames defined at
! the beginning of the forcing interval.
! filename_at_begin_of_interval = .false. for filenames defined at
! the middle of the forcing interval.
!
!-----------------------------------------------------------------------
filename_at_begin_of_interval = .false.
!-----------------------------------------------------------------------
!
! find the time associated with the beginning or middle of the
! forcing interval.
!
!-----------------------------------------------------------------------
if (filename_at_begin_of_interval) then
file_time = forcing_time - p5*forcing_data_inc + eps
else ! filename defined at forcing time (middle of forcing interval)
file_time = forcing_time + eps
endif
!-----------------------------------------------------------------------
!
! if the file time is negative, find the first appropriate
! positive value and use it.
!
!-----------------------------------------------------------------------
if (file_time < c0 ) then
if (my_task == master_task) then
write(stdout,'(a52)') &
'WARNING: apparently need a forcing file earlier than'
write(stdout,'(a54)') &
'01-01-0000 will attempt to find first appropriate file'
endif
do n = 1,10
if (filename_at_begin_of_interval) then
file_time = forcing_time + (n-p5)*forcing_data_inc + eps
else
file_time = forcing_time + n*forcing_data_inc + eps
endif
if(file_time >= c0 ) exit
enddo
endif
if (file_time < c0 ) then
call exit_POP(sigAbort,'get_forcing_filename confused')
endif
!-----------------------------------------------------------------------
!
! convert file time to time relative to the beginning of the year.
! then find the day and hour corresponding to that time.
!
! note that the use of 'eps' above is important here due to the
! 'int' and 'nint' functions, especially if file_time00 < 0.
!
!-----------------------------------------------------------------------
thour00_begin_this_year = thour00 - &
(elapsed_days_this_year + frac_day)*24.0_r8
file_time00 = file_time - thour00_begin_this_year
file_day = int(abs(file_time00)/24.0_r8)
file_hour = nint(abs(file_time00) - 24*file_day)
!-----------------------------------------------------------------------
!
! if file_time00 < 0, then the time is in the previous year.
!
! if file_time00 > the number of hours in a year, then the time is
! in the next year.
!
!-----------------------------------------------------------------------
if (file_time00 < c0 )then
file_hour = 24 - file_hour ! counting backwards
file_year = iyear - 1
days_in_prev_year = days_in_norm_year
if (allow_leapyear .and. mod(file_year,4) == 0) then ! check if previous
if (mod(file_year,100) /= 0) & ! year was leapyear
days_in_prev_year = days_in_leap_year
if (mod(file_year,400) == 0) &
days_in_prev_year = days_in_leap_year
endif
file_day = days_in_prev_year - file_day
elseif (file_time00 > 24*days_in_year) then
file_day = file_day - days_in_year + 1
file_year = iyear + 1
else
file_day = file_day + 1 ! day counting starts at 1, not 0
file_year = iyear
endif
!-----------------------------------------------------------------------
!
! write year, day, and hour into character variables and use them
! to create filename.
!
!-----------------------------------------------------------------------
write(file_chour,'(i3)') 100 + file_hour
write(file_cday, '(i4)') 1000 + file_day
write(file_cyear,'(i5)') 10000 + file_year
forcing_filename = char_blank
cindx = len_trim(forcing_infile)
forcing_filename(1:cindx) = forcing_infile(1:cindx)
cindx = cindx + 1
forcing_filename(cindx:cindx) = '.'
cindx = cindx + 1
forcing_filename(cindx:cindx+3) = file_cyear(2:5)
cindx = cindx + 4
forcing_filename(cindx:cindx) = '.'
cindx = cindx + 1
forcing_filename(cindx:cindx+2) = file_cday(2:4)
cindx = cindx + 3
forcing_filename(cindx:cindx) = '.'
cindx = cindx + 1
forcing_filename(cindx:cindx+1) = file_chour(2:3)
!-----------------------------------------------------------------------
!EOC
end subroutine get_forcing_filename
!***********************************************************************
!BOP
! !IROUTINE: find_forcing_times
! !INTERFACE:
subroutine find_forcing_times(forcing_time, forcing_data_inc, & 25,1
forcing_interp_type, forcing_time_next, &
forcing_time_min_loc, forcing_data_update, &
forcing_data_type)
! !DESCRIPTION:
! Find the times associated with forcing data needed for
! interpolation given the current time and time increment
! between forcing data fields. forcing times are assumed
! to be centered; for example, with daily forcing, the value
! is assumed to occur at noon (hour = 12).
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character (char_len), intent(in) :: &
forcing_data_type, &! frequency of forcing data
forcing_interp_type ! type of interpolation to be used
real (r8), intent(in) :: &
forcing_data_inc ! time increment between forcing data
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(12), intent(inout) :: &
forcing_time
! !OUTPUT PARAMETERS:
integer (int_kind), intent(out) :: &
forcing_time_min_loc ! index to first location (in time)
! of forcing data array
real (r8), intent(out) :: &
forcing_time_next, &! next time for updating forcing
forcing_data_update ! next time to update forcing data
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
nstart, n, &
second, third, fourth
real (r8) :: &
forcing_time_test, &
forcing_time_min
!-----------------------------------------------------------------------
!
! set mid_month and end_month times for monthly options.
!
!-----------------------------------------------------------------------
if (forcing_data_type == 'monthly-calendar') then
thour00_midmonth = thour00_midmonth_calendar
thour00_endmonth = thour00_endmonth_calendar
else if (forcing_data_type == 'monthly-equal') then
thour00_midmonth = thour00_midmonth_equal
thour00_endmonth = thour00_endmonth_equal
endif
!-----------------------------------------------------------------------
!
! for 'monthly-equal' and 'monthly-calendar' set the forcing times to
! be the mid-month times plus the time at the beginning of the year.
! then find the first forcing time greater than the current time
! and exit. note that the loop will not exit if the current time
! is in the second half of December in which case n will have the
! value 13.
! for 'n-hour' forcing, start at a time that is definitely less
! than the first forcing time greater than the current time and
! just increment the test time by the forcing increment until we
! find a forcing time greater than the current time.
!
!-----------------------------------------------------------------------
select case (forcing_data_type)
case ('monthly-equal','monthly-calendar')
forcing_time(1:12) = thour00_begin_this_year + &
thour00_midmonth(1:12)
do n = 1,12
if (forcing_time(n) > thour00) exit
enddo
case ('n-hour')
nstart = int(thour00/forcing_data_inc) - 3
do n = nstart, nstart + 6
forcing_time_test = (n - p5)*forcing_data_inc
if (forcing_time_test > thour00) exit
enddo
end select
!-----------------------------------------------------------------------
!
! Set forcing times depending on the type of temporal interpolation.
! Also set forcing_time_min_loc, which is the temporal index of the
! minimum forcing time to be used for the interpolation on the forcing
! data arrays. for 'monthly' data, forcing_time_min_loc varies between
! 1 and 12. For example, if a 4point interpolation is required and
! the months that are needed are 3,4,5,6 then forcing_time_min_loc = 3.
! For 'n-hour' data, forcing_time_min_loc = 1 always because
! initially the data is read chronologically 1 file at a time.
! Also set the times for the next interpolation and next time
! that new data will be needed.
!
!-----------------------------------------------------------------------
select case(forcing_interp_type)
case ('nearest') ! select nearest forcing time
select case (forcing_data_type)
case ('monthly-equal','monthly-calendar')
do n = 1,12
if (thour00 <= &
thour00_begin_this_year + thour00_endmonth(n)) exit
enddo
if (n == 13) call exit_POP(sigAbort, &
'Error finding current month in find_forcing_times')
forcing_time_min_loc = n
forcing_time_min_loc = mod(forcing_time_min_loc,12)
if (forcing_time_min_loc <= 0) &
forcing_time_min_loc = forcing_time_min_loc + 12
if (n == 12) then
forcing_time_next = forcing_time(1) + hours_in_year ! next Jan
else
forcing_time_next = forcing_time(forcing_time_min_loc + 1)
endif
forcing_data_update = thour00_begin_this_year + &
thour00_endmonth(forcing_time_min_loc)
case ('n-hour')
if ((forcing_time_test-thour00) <= p5*forcing_data_inc ) then
forcing_time_min = (n - p5)*forcing_data_inc
else
forcing_time_min = (n - 1.5_r8)*forcing_data_inc
endif
forcing_time_min_loc = 1
forcing_time(forcing_time_min_loc) = forcing_time_min
forcing_time_next = forcing_time(forcing_time_min_loc) + &
forcing_data_inc
forcing_data_update = forcing_time(forcing_time_min_loc) + &
p5*forcing_data_inc
end select
case ('linear') ! linear interpolation
select case(forcing_data_type)
case ('monthly-equal','monthly-calendar')
forcing_time_min_loc = n - 1
forcing_time_min_loc = mod(forcing_time_min_loc,12)
if (forcing_time_min_loc <= 0) &
forcing_time_min_loc = forcing_time_min_loc + 12
second = mod(forcing_time_min_loc+1,12)
if (second == 0) second = 12
if (n == 1) then
forcing_time(forcing_time_min_loc) = &
forcing_time(forcing_time_min_loc) - hours_in_year ! previous Dec
elseif (n >= 3) then
forcing_time(1:n-2) = forcing_time(1:n-2) + hours_in_year
endif
if (n == 12) then
forcing_time_next = forcing_time(1) ! next Jan
else
forcing_time_next = forcing_time(second+1) ! next month
endif
case ('n-hour')
forcing_time_min_loc = 1
second = 2
forcing_time(forcing_time_min_loc) = (n - 1.5_r8)* &
forcing_data_inc
forcing_time(second) = forcing_time(forcing_time_min_loc) + &
forcing_data_inc
forcing_time_next = forcing_time(forcing_time_min_loc) + &
c2*forcing_data_inc
end select
forcing_data_update = forcing_time(second)
case ('4point') ! 4-point interpolation
select case(forcing_data_type)
case ('monthly-equal','monthly-calendar')
forcing_time_min_loc = n - 2
forcing_time_min_loc = mod(forcing_time_min_loc,12)
if (forcing_time_min_loc <= 0) &
forcing_time_min_loc = forcing_time_min_loc + 12
second = mod(forcing_time_min_loc+1,12)
third = mod(forcing_time_min_loc+2,12)
fourth = mod(forcing_time_min_loc+3,12)
if (second == 0) second = 12
if (third == 0) third = 12
if (fourth == 0) fourth = 12
if (n == 1) then
forcing_time(forcing_time_min_loc) = &
forcing_time(forcing_time_min_loc) - hours_in_year ! previous Nov
forcing_time(second) = &
forcing_time(second) - hours_in_year ! previous Dec
else if (n == 2) then
forcing_time(forcing_time_min_loc) = &
forcing_time(forcing_time_min_loc) - hours_in_year ! previous Dec
else if (n >= 4) then
forcing_time(1:n-3) = forcing_time(1:n-3) + hours_in_year
endif
if (n == 11) then
forcing_time_next = forcing_time(1) ! next Jan
else
forcing_time_next = forcing_time(fourth+1) ! next month
endif
case ('n-hour')
forcing_time_min_loc = 1
second = 2
third = 3
fourth = 4
forcing_time(forcing_time_min_loc) = (n - 2.5_r8)* &
forcing_data_inc
forcing_time(second) = forcing_time(forcing_time_min_loc) + &
forcing_data_inc
forcing_time(third) = forcing_time(forcing_time_min_loc) + &
c2*forcing_data_inc
forcing_time(fourth) = forcing_time(forcing_time_min_loc) + &
c3*forcing_data_inc
forcing_time_next = forcing_time(forcing_time_min_loc) + &
c4*forcing_data_inc
end select
forcing_data_update = forcing_time(third)
end select
!-----------------------------------------------------------------------
!EOC
end subroutine find_forcing_times
!***********************************************************************
!BOC
! !IROUTINE: interpolate_forcing
! !INTERFACE:
subroutine interpolate_forcing(INTERP, FIELD, forcing_time, & 25,3
forcing_interp_type, forcing_time_min_loc, &
forcing_interp_freq, forcing_interp_inc, forcing_interp_next, &
forcing_interp_last, nsteps_run_check)
! !DESCRIPTION:
! Temporally interpolates field based on input interpolation options
! and indices.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character (char_len), intent(in) :: &
forcing_interp_type, &! type of interpolation to use
forcing_interp_freq ! frequency for interpolating forcing
integer (int_kind), intent(in) :: &
nsteps_run_check, &! check for handling first time step
forcing_time_min_loc ! first location of data to interpolate
real (r8), dimension(:,:,:,:,:), intent(in) :: &
FIELD ! forcing field to be interpolated
real (r8), intent(in) :: &
forcing_interp_inc, &! forcing increment
forcing_interp_next
real (r8), dimension(12), intent(in) :: &
forcing_time ! times at which forcing located for interpolation
! !INPUT/OUTPUT PARAMETERS:
real (r8), intent(inout) :: &
forcing_interp_last ! time when last interpolation done
! !OUTPUT PARAMETERS:
real (r8), dimension(:,:,:,:), intent(out) :: &
INTERP ! result of interpolation
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k,n, &! loop indices
iblock, &! dummy block index
second, third, &! time indices for data to be interpolated
fourth, &!
ntime ! number of time slices in data array
real (r8) :: &
weight ! linear interpolation weight
!-----------------------------------------------------------------------
!
! if interpolation is done every timestep, then interpolate to
! the current time.
! if interpolation is done every n-hours, then interpolate to
! the current interpolation time (which is >= current time)
! or if it is the first step of the current sub-run, the value
! of forcing_interp_last was read in from the restart file.
!
!-----------------------------------------------------------------------
ntime = size(FIELD, dim=5)
if (forcing_interp_freq == 'every-timestep') then
forcing_interp_last = thour00
else if (forcing_interp_freq == 'n-hour') then
if (nsteps_run_check /= 0) forcing_interp_last = &
forcing_interp_next
else
call exit_POP(sigAbort,'Error in interpolate_forcing')
endif
select case(forcing_interp_type)
case ('nearest')
INTERP = FIELD(:,:,:,:,forcing_time_min_loc)
case ('linear')
second = mod(forcing_time_min_loc+1,ntime)
if (second == 0) second = ntime
weight = (forcing_time(second) - forcing_interp_last)/ &
(forcing_time(second) - &
forcing_time(forcing_time_min_loc))
!*** 3d interior forcing and surface forcing have different
!*** order of indices so must be treated differently
if (size(FIELD,dim=3) == km) then ! 3d interior forcing
!$OMP PARALLEL DO
do iblock=1,nblocks_clinic
INTERP(:,:,:,iblock) = &
weight *FIELD(:,:,:,iblock,forcing_time_min_loc) + &
(c1 - weight)*FIELD(:,:,:,iblock,second)
end do
!$OMP END PARALLEL DO
else ! 2d surface forcing with 1 or more fields
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,size(FIELD,dim=4)
INTERP(:,:,iblock,n) = &
weight *FIELD(:,:,iblock,n,forcing_time_min_loc) + &
(c1 - weight)*FIELD(:,:,iblock,n,second)
end do
end do
!$OMP END PARALLEL DO
endif
case ('4point')
second = mod(forcing_time_min_loc+1,ntime)
third = mod(forcing_time_min_loc+2,ntime)
fourth = mod(forcing_time_min_loc+3,ntime)
if (second == 0) second = ntime
if (third == 0) third = ntime
if (fourth == 0) fourth = ntime
!*** 3d interior forcing and surface forcing have different
!*** order of indices so must be treated differently
if (size(FIELD,dim=3) == km) then ! 3d interior forcing
!$OMP PARALLEL DO PRIVATE(iblock, k)
do iblock=1,nblocks_clinic
do k=1,km
call interp_4pt(INTERP(:,:,k,iblock), &
FIELD(:,:,k,iblock,forcing_time_min_loc), &
FIELD(:,:,k,iblock,second), &
FIELD(:,:,k,iblock,third), &
FIELD(:,:,k,iblock,fourth), &
forcing_time(forcing_time_min_loc), &
forcing_time(second),forcing_time(third), &
forcing_time(fourth),forcing_interp_last)
end do
end do
!$OMP END PARALLEL DO
else ! 2d surface forcing with 1 or more fields
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,size(FIELD,dim=4)
call interp_4pt(INTERP(:,:,iblock,n), &
FIELD(:,:,iblock,n,forcing_time_min_loc), &
FIELD(:,:,iblock,n,second), &
FIELD(:,:,iblock,n,third), &
FIELD(:,:,iblock,n,fourth), &
forcing_time(forcing_time_min_loc), &
forcing_time(second),forcing_time(third), &
forcing_time(fourth),forcing_interp_last)
end do
end do
!$OMP END PARALLEL DO
endif
end select
!-----------------------------------------------------------------------
!EOC
end subroutine interpolate_forcing
!***********************************************************************
!BOP
! !IROUTINE:
! !INTERFACE:
subroutine update_forcing_data(forcing_time, forcing_time_min_loc, & 25,16
forcing_interp_type, forcing_data_next, &
forcing_data_update, forcing_data_type, &
forcing_data_inc, FIELD, &
forcing_data_rescale, &
forcing_data_label, forcing_data_names, &
forcing_bndy_loc, forcing_bndy_type, &
forcing_infile, forcing_infile_fmt)
! !DESCRIPTION:
! Updates the data to be used in the next temporal interpolation.
! If data is monthly this is a shuffling of indices since
! all 12 months are in memory. if data is n-hour, new data
! needs to be read in as well as shuffling of indices.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character (*), intent(in) :: &
forcing_infile, &! file containing forcing data
forcing_infile_fmt, &! format (bin or netcdf) for input file
forcing_data_type, &! type of forcing data
forcing_interp_type, &! interpolation to use for forcing
forcing_data_label ! description of forcing
character (char_len), dimension(:), intent(in) :: &
forcing_data_names ! short names for input forcing fields
integer (int_kind), dimension(:), intent(in) :: &
forcing_bndy_loc, &! location and field type for ghost
forcing_bndy_type ! cell updates
real (r8), dimension(20), intent(in) :: &
forcing_data_rescale ! factors for converting to model units
real (r8), intent(in) :: &
forcing_data_inc ! time increment between forcing data
! !INPUT/OUTPUT PARAMETERS:
integer (int_kind), intent(inout) :: &
forcing_time_min_loc ! first temporal index for interpolation
real (r8), dimension(12), intent(inout) :: &
forcing_time !
real (r8), intent(inout) :: &
forcing_data_next, &! time next data required
forcing_data_update ! time data to be updated
real (r8), dimension(:,:,:,:,:), intent(inout) :: &
FIELD ! forcing field to be read if necessary
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
n, &! dummy index
forcing_time_max_loc ! time index for last forcing time
character (char_len) :: &
forcing_filename ! name of forcing file
type (datafile) :: &
in_forcing_file ! data file type for input forcing file
type (io_field_desc) :: &
forcing_field ! io type for input forcing field
type (io_dim) :: &
i_dim, j_dim, &! dimension descriptors for horiz dims
k_dim ! dimension descriptor for vertical levels
!-----------------------------------------------------------------------
!
! determine data by forcing type
!
!-----------------------------------------------------------------------
select case(forcing_data_type)
!-----------------------------------------------------------------------
!
! for monthly forcing, all twelve months are stored - simply compute
! proper indices into forcing data array
!
!-----------------------------------------------------------------------
case ('monthly-equal','monthly-calendar')
!*** compute first forcing time
forcing_time(forcing_time_min_loc) = &
forcing_time(forcing_time_min_loc) + hours_in_year
forcing_time_min_loc = forcing_time_min_loc + 1
if (forcing_time_min_loc > 12) forcing_time_min_loc = 1
!*** compute times for other times for use in interpolation
select case(forcing_interp_type)
case ('nearest')
forcing_time_max_loc = forcing_time_min_loc
forcing_data_update = forcing_data_update + &
24.0_r8*days_in_month(forcing_time_max_loc)
case ('linear')
forcing_time_max_loc = mod(forcing_time_min_loc + 1,12)
if (forcing_time_max_loc == 0) forcing_time_max_loc = 12
forcing_data_update = forcing_time(forcing_time_max_loc)
case ('4point')
forcing_time_max_loc = mod(forcing_time_min_loc + 3,12)
if (forcing_time_max_loc == 0) forcing_time_max_loc = 12
n = forcing_time_max_loc - 1
if (n == 0) n = 12
forcing_data_update = forcing_time(n)
end select
!*** compute next forcing time
forcing_time(forcing_time_max_loc) = forcing_data_next
n = mod(forcing_time_max_loc + 1, 12)
if (n == 0) n = 12
forcing_data_next = forcing_time(n)
!-----------------------------------------------------------------------
!
! for n-hour option, read in required forcing data from forcing file
!
!-----------------------------------------------------------------------
case ('n-hour')
call get_forcing_filename(forcing_filename, forcing_infile, &
forcing_data_next, forcing_data_inc)
in_forcing_file = construct_file(forcing_infile_fmt, &
full_name=trim(forcing_filename), &
record_length=rec_type_dbl, &
recl_words=nx_global*ny_global)
call data_set(in_forcing_file, 'open_read')
i_dim = construct_io_dim('i',nx_global)
j_dim = construct_io_dim('j',ny_global)
if (size(FIELD,dim=3) == km) then ! 3d interior forcing field
k_dim = construct_io_dim('k',km)
forcing_field = construct_io_field( &
trim(forcing_data_names(1)), &
dim1=i_dim, dim2=j_dim, dim3=k_dim, &
field_loc = forcing_bndy_loc(1), &
field_type = forcing_bndy_type(1), &
d3d_array=FIELD(:,:,:,:,forcing_time_min_loc))
call data_set (in_forcing_file, 'define', forcing_field)
call data_set (in_forcing_file, 'read', forcing_field)
call destroy_io_field(forcing_field)
!*** renormalize if necessary to compensate for different
!*** units.
if (forcing_data_rescale(1) /= c1) &
FIELD(:,:,:,:,forcing_time_min_loc) = &
FIELD(:,:,:,:,forcing_time_min_loc)*forcing_data_rescale(1)
else ! forcing field(s) for surface forcing
do n = 1,size(FIELD,dim=4)
forcing_field = construct_io_field( &
trim(forcing_data_names(n)), &
dim1=i_dim, dim2=j_dim, &
field_loc = forcing_bndy_loc(n), &
field_type = forcing_bndy_type(n), &
d2d_array=FIELD(:,:,:,n,forcing_time_min_loc))
call data_set (in_forcing_file, 'define', forcing_field)
call data_set (in_forcing_file, 'read', forcing_field)
call destroy_io_field(forcing_field)
!*** renormalize if necessary to compensate for different
!*** units.
if (forcing_data_rescale(n) /= c1) &
FIELD(:,:,:,n,forcing_time_min_loc) = &
FIELD(:,:,:,n,forcing_time_min_loc)*forcing_data_rescale(n)
enddo
endif
call data_set (in_forcing_file, 'close')
call destroy_file(in_forcing_file)
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,'(a9,f12.3,1x,a,a12,a)') 'tday00 = ', tday00, &
trim(forcing_data_label),' file read: ', &
trim(forcing_filename)
endif
forcing_time(forcing_time_min_loc) = forcing_data_next
forcing_time_min_loc = forcing_time_min_loc + 1
if (forcing_time_min_loc > size(FIELD,dim=5)) &
forcing_time_min_loc = 1
forcing_data_next = forcing_data_next + forcing_data_inc
forcing_data_update = forcing_data_update + forcing_data_inc
!-----------------------------------------------------------------------
end select
!-----------------------------------------------------------------------
!EOC
end subroutine update_forcing_data
!***********************************************************************
!BOP
! !IROUTINE: echo_forcing_options
! !INTERFACE:
subroutine echo_forcing_options(forcing_data_type, & 6
forcing_formulation, &
forcing_data_inc, &
forcing_interp_freq, &
forcing_interp_type, &
forcing_interp_inc, &
forcing_label)
! !DESCRIPTION:
! Writes out forcing options to stdout (or log file).
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), intent(in) :: &
forcing_data_inc, &
forcing_interp_inc
character (char_len), intent(in) :: &
forcing_data_type, &
forcing_interp_freq, &
forcing_interp_type, &
forcing_label, &
forcing_formulation
!EOP
!BOC
!-----------------------------------------------------------------------
!
! echo forcing selections to stdout. only write from master task.
!
!-----------------------------------------------------------------------
if (my_task == master_task) then
!-----------------------------------------------------------------------
!
! no forcing
!
!-----------------------------------------------------------------------
if (forcing_data_type == 'none') then
write(stdout,'(a3,a)') 'No ', trim(forcing_label)
!-----------------------------------------------------------------------
!
! if there is some kind of forcing, write out the type (eg, annual)
! and the formulation (eg, restoring).
! then write increment between forcing data if appropriate.
!
!-----------------------------------------------------------------------
else
if (len_trim(forcing_formulation) == 0) then
! Work around current Linux runtime error
write(stdout,'(a,a1,a)') trim(forcing_label), ':', &
trim(forcing_data_type)
else
write(stdout,'(a,a1,a,a1,a)') trim(forcing_label), ':', &
trim(forcing_data_type),' ', &
trim(forcing_formulation)
end if
if (forcing_data_type == 'n-hour') then
write(stdout,'(a,a29,f7.3)') trim(forcing_label), &
' Data Time Increment (days): ', &
forcing_data_inc/24.0_r8
endif
!-----------------------------------------------------------------------
!
! write the frequency that the forcing is updated.
!
!-----------------------------------------------------------------------
select case(forcing_interp_freq)
case ('never')
write(stdout,'(a,a27)') trim(forcing_label), &
' data is never interpolated'
case ('n-hour')
write(stdout,'(a,a25,f7.3,a5)') trim(forcing_label), &
' data interpolated every ', &
forcing_interp_inc/24.0_r8,' days'
case ('every-timestep')
write(stdout,'(a,a33)') trim(forcing_label), &
' data interpolated every timestep'
end select
!-----------------------------------------------------------------------
!
! write order of interpolation if appropriate.
!
!-----------------------------------------------------------------------
if (forcing_data_type /= 'none' .and. &
forcing_data_type /= 'analytic' .and. &
forcing_data_type /= 'annual') &
write(stdout,'(a,a21,a)') trim(forcing_label), &
' Interpolation type: ',trim(forcing_interp_type)
endif ! forcing type /= none
write(stdout,blank_fmt)
!-----------------------------------------------------------------------
!
! end of output
!
!-----------------------------------------------------------------------
endif ! master task
!-----------------------------------------------------------------------
!EOC
end subroutine echo_forcing_options
!***********************************************************************
!BOP
! !IROUTINE: interp_4pt
! !INTERFACE:
subroutine interp_4pt(INTERP, DATA1, DATA2, DATA3, DATA4, & 2,6
tdata1, tdata2, tdata3, tdata4, time)
! !DESCRIPTION:
! Interpolates data at four consecutive (in time) points to a
! particular time using a 3rd degree polynomial fit. This is an
! implementation of Nevilles algorithm as described in Numerical
! Recipes.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(in) :: &
DATA1, DATA2, DATA3, DATA4 ! data for each grid point in
! a block for 4 consecutive
! instants in time
real (r8), intent(in) :: &
time, &! time for which interpolated data required
tdata1, tdata2, &! time points at which above data is located
tdata3, tdata4
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
INTERP ! data interpolated to input time
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
real (r8), dimension(nx_block,ny_block) :: &
WORK1,WORK2,WORK3,WORK4 ! local work space
!-----------------------------------------------------------------------
WORK1 = c0
WORK2 = c0
WORK3 = c0
WORK4 = c0
call det(WORK1,DATA1,DATA2,tdata1,tdata2,time) ! P12
call det(WORK2,DATA2,DATA3,tdata2,tdata3,time) ! P23
call det(WORK3,DATA3,DATA4,tdata3,tdata4,time) ! P34
call det(WORK4,WORK1,WORK2,tdata1,tdata3,time) ! P123
call det(WORK1,WORK2,WORK3,tdata2,tdata4,time) ! P234
call det(INTERP,WORK4,WORK1,tdata1,tdata4,time) ! P1234
!-----------------------------------------------------------------------
!EOC
end subroutine interp_4pt
!***********************************************************************
!BOP
! !IROUTINE: det
! !INTERFACE:
subroutine det(C,A,B,y,z,x) 6
! !DESCRIPTION:
! Determinant routine used by Nevilles algorithm in the 4-point
! time interpolation.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), intent(in) :: &
y,z,x
real (r8), dimension(nx_block,ny_block), intent(in) :: &
A,B
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
C ! output determinant
!-----------------------------------------------------------------------
C = ( A*(z-x) - B*(y-x) )/(z-y)
!-----------------------------------------------------------------------
!EOC
end subroutine det
!***********************************************************************
end module forcing_tools
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||