diff --git a/.gitmodules b/.gitmodules index 692d9223..949b88cd 100644 --- a/.gitmodules +++ b/.gitmodules @@ -14,7 +14,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/ESCOMP/atmospheric_physics - fxtag = 098585940ad763be58ebab849bb8eaf325fda42a + fxtag = atmos_phys0_05_000 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index c12d23a1..ba8fbe40 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -333,6 +333,19 @@ + + + char*256 + tropopause + tropopause_nl + + Full pathname of boundary dataset for tropopause climatology. + + + ${DIN_LOC_ROOT}/atm/cam/chem/trop_mozart/ub/clim_p_trop.nc + + + @@ -364,7 +377,7 @@ cam_logfile_nl 0,1,2,3 - Include exra checks and logging output. + Include extra checks and logging output. 0: No extra checks or output 1: Extra informational output 2: Level 1 plus extra checks and informational output diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 0782738a..0f432c5e 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -23,10 +23,12 @@ module cam_comp use time_manager, only: timemgr_init, get_step_size use time_manager, only: get_nstep use time_manager, only: is_first_step, is_first_restart_step + use time_manager, only: get_curr_calday use camsrfexch, only: cam_out_t, cam_in_t use physics_types, only: phys_state, phys_tend use physics_types, only: dtime_phys + use physics_types, only: calday use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -98,6 +100,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & use physics_grid, only: columns_on_task use vert_coord, only: pver use phys_vars_init_check, only: mark_as_initialized + use tropopause_climo_read, only: tropopause_climo_read_file ! Arguments character(len=cl), intent(in) :: caseid ! case ID @@ -163,11 +166,16 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & call cam_ctrl_set_orbit(eccen, obliqr, lambm0, mvelpp) + call timemgr_init( & dtime, calendar, start_ymd, start_tod, ref_ymd, & ref_tod, stop_ymd, stop_tod, curr_ymd, curr_tod, & perpetual_run, perpetual_ymd, initial_run_in) + ! Get current fractional calendar day. Needs to be updated at every timestep. + calday = get_curr_calday() + call mark_as_initialized('fractional_calendar_days_on_end_of_current_timestep') + ! Read CAM namelists. filein = "atm_in" // trim(inst_suffix) call read_namelist(filein, single_column, scmlat, scmlon) @@ -224,6 +232,9 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & !!XXgoldyXX: ^ need to import this end if + ! Read tropopause climatology + call tropopause_climo_read_file() + call phys_init() !!XXgoldyXX: v need to import this @@ -270,6 +281,9 @@ subroutine cam_timestep_init() ! call phys_timestep_init() + ! Update current fractional calendar day. Needs to be updated at every timestep. + calday = get_curr_calday() + end subroutine cam_timestep_init ! !----------------------------------------------------------------------- diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 12081a81..7031be1f 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -44,6 +44,8 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) ! use cam_diagnostics, only: diag_readnl use inic_analytic_utils, only: analytic_ic_readnl + use tropopause_climo_read, only: tropopause_climo_readnl + ! use tracers, only: tracers_readnl ! use nudging, only: nudging_readnl @@ -99,6 +101,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) ! call diag_readnl(nlfilename) ! call check_energy_readnl(nlfilename) call analytic_ic_readnl(nlfilename) + call tropopause_climo_readnl(nlfilename) ! call scam_readnl(nlfilename, single_column, scmlat, scmlon) ! call nudging_readnl(nlfilename) diff --git a/src/data/generate_registry_data.py b/src/data/generate_registry_data.py index 33adc954..187cf7be 100755 --- a/src/data/generate_registry_data.py +++ b/src/data/generate_registry_data.py @@ -26,6 +26,7 @@ sys.path.append(__SPINSCRIPTS) # end if _ALL_STRINGS_REGEX = re.compile(r'[A-Za-z][A-Za-z_0-9]+') +_FORTRAN_NUMERIC_REGEX = re.compile(r'^[+-]?(\d+\.?\d*|\.\d+)([eE][+-]?\d+)?(_kind_phys)?$') # CCPP framework imports # pylint: disable=wrong-import-position @@ -940,9 +941,9 @@ def add_variable(self, newvar): all_strings = _ALL_STRINGS_REGEX.findall(newvar.initial_value) init_val_vars = set() excluded_initializations = {'null', 'nan', 'false', 'true'} - # Exclude NULL and nan variables + # Exclude NULL and nan variables and valid Fortran numeric values that pass the string regex (e.g. 1.e36, -3.2e5) for var in all_strings: - if var.lower() not in excluded_initializations: + if var.lower() not in excluded_initializations and not _FORTRAN_NUMERIC_REGEX.match(newvar.initial_value): init_val_vars.add(var) # end if # end if diff --git a/src/data/registry.xml b/src/data/registry.xml index 9bca0f62..ef2abbc2 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -10,6 +10,7 @@ $SRCROOT/src/data/physconst.meta $SRCROOT/src/physics/utils/physics_grid.meta $SRCROOT/src/physics/utils/cam_constituents.meta + $SRCROOT/src/physics/utils/tropopause_climo_read.meta $SRCROOT/src/data/air_composition.meta $SRCROOT/src/data/cam_thermo.meta $SRCROOT/src/data/ref_pres.meta @@ -348,6 +349,12 @@ current timestep number 0 + + + fractional calendar day at end of current timestep + ccpp error message '' + + + Fill value for diagnostic/history outputs + 1.e36_kind_phys + \section arg_table_tropopause_climo_read Argument Table +!! \htmlinclude tropopause_climo_read.html + ! months in year for climatological tropopause pressure data + integer, public, parameter :: tropp_slices = 12 + + ! climatological tropopause pressures (ncol,ntimes) + real(kind_phys), public, allocatable :: tropp_p_loc(:,:) + + ! monthly day-of-year times corresponding to climatological data (12) + real(kind_phys), public, allocatable :: tropp_days(:) + + ! Private module data + character(len=shr_kind_cl) :: tropopause_climo_file = unset_str + +contains + ! Read namelist variable tropopause_climo_file. + ! Containing this within CAM-SIMA instead of within scheme as otherwise the climo filepath + ! has to be passed from physics -> here then the data from here -> physics. + subroutine tropopause_climo_readnl(nlfile) + use shr_nl_mod, only: find_group_name => shr_nl_find_group_name + use shr_kind_mod, only: shr_kind_cm + use mpi, only: mpi_character + use spmd_utils, only: mpicom + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use spmd_utils, only: masterproc + + ! filepath for file containing namelist input + character(len=*), intent(in) :: nlfile + + ! Local variables + integer :: unitn, ierr + character(len=*), parameter :: subname = 'tropopause_climo_readnl' + character(len=shr_kind_cm) :: errmsg + + namelist /tropopause_nl/ tropopause_climo_file + + errmsg = '' + + if (masterproc) then + open(newunit=unitn, file=trim(nlfile), status='old') + call find_group_name(unitn, 'tropopause_nl', status=ierr) + if (ierr == 0) then + read(unitn, tropopause_nl, iostat=ierr, iomsg=errmsg) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading namelist:' // errmsg) + end if + end if + close(unitn) + end if + + ! Broadcast namelist variables + call mpi_bcast(tropopause_climo_file, len(tropopause_climo_file), mpi_character, 0, mpicom, ierr) + + ! Print out namelist variables + if (masterproc) then + write(iulog,*) subname, ' options:' + write(iulog,*) ' Tropopause climatology file will be read: ', trim(tropopause_climo_file) + endif + + end subroutine tropopause_climo_readnl + + subroutine tropopause_climo_read_file() + !------------------------------------------------------------------ + ! ... read tropopause climatology dataset file + !------------------------------------------------------------------ + use shr_kind_mod, only: shr_kind_cm + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use spmd_utils, only: masterproc + use interpolate_data, only: lininterp_init, lininterp, interp_type, lininterp_finish + use physics_grid, only: get_rlat_all_p, get_rlon_all_p + use physics_grid, only: pcols => columns_on_task + use physconst, only: pi + use time_manager, only: get_calday + use ioFileMod, only: cam_get_file + use cam_pio_utils, only: cam_pio_openfile + use pio, only: file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen + use pio, only: pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite + use phys_vars_init_check, only: mark_as_initialized + + !------------------------------------------------------------------ + ! ... local variables + !------------------------------------------------------------------ + integer :: i, j, n + integer :: ierr + type(file_desc_t) :: pio_id + integer :: dimid + type(var_desc_t) :: vid + integer :: nlon, nlat, ntimes + integer :: start(3) + integer :: count(3) + + ! YMD (01-16, 02-14, ...) corresponding to the time slices of the climatological + ! tropopause dataset. Will be converted to day-of-year and stored in tropp_days. + integer, parameter :: dates(tropp_slices) = (/ 116, 214, 316, 415, 516, 615, & + 716, 816, 915, 1016, 1115, 1216 /) + type(interp_type) :: lon_wgts, lat_wgts + real(kind_phys), allocatable :: tropp_p_in(:,:,:) + real(kind_phys), allocatable :: lat(:) + real(kind_phys), allocatable :: lon(:) + real(kind_phys) :: to_lats(pcols), to_lons(pcols) + real(kind_phys), parameter :: d2r=pi/180._kind_phys, zero=0._kind_phys, twopi=pi*2._kind_phys + character(len=shr_kind_cl) :: locfn + character(len=shr_kind_cm) :: errmsg + + errmsg = '' + + !----------------------------------------------------------------------- + ! ... open netcdf file + !----------------------------------------------------------------------- + call cam_get_file (tropopause_climo_file, locfn, allow_fail=.false.) + call cam_pio_openfile(pio_id, trim(locfn), PIO_NOWRITE) + + !----------------------------------------------------------------------- + ! ... get time dimension + !----------------------------------------------------------------------- + ierr = pio_inq_dimid( pio_id, 'time', dimid ) + ierr = pio_inq_dimlen( pio_id, dimid, ntimes ) + if( ntimes /= tropp_slices ) then + write(iulog,*) 'tropopause_climo_read_file: number of months = ',ntimes,'; expecting ',tropp_slices + call endrun('tropopause_climo_read_file: number of months in file incorrect') + end if + !----------------------------------------------------------------------- + ! ... get latitudes + !----------------------------------------------------------------------- + ierr = pio_inq_dimid( pio_id, 'lat', dimid ) + ierr = pio_inq_dimlen( pio_id, dimid, nlat ) + allocate( lat(nlat), stat=ierr, errmsg=errmsg ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_climo_read_file: lat allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate lat, error = ' // errmsg) + end if + ierr = pio_inq_varid( pio_id, 'lat', vid ) + ierr = pio_get_var( pio_id, vid, lat ) + lat(:nlat) = lat(:nlat) * d2r + !----------------------------------------------------------------------- + ! ... get longitudes + !----------------------------------------------------------------------- + ierr = pio_inq_dimid( pio_id, 'lon', dimid ) + ierr = pio_inq_dimlen( pio_id, dimid, nlon ) + allocate( lon(nlon), stat=ierr, errmsg=errmsg ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_climo_read_file: lon allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate lon, error = ' // errmsg) + end if + ierr = pio_inq_varid( pio_id, 'lon', vid ) + ierr = pio_get_var( pio_id, vid, lon ) + lon(:nlon) = lon(:nlon) * d2r + + !------------------------------------------------------------------ + ! ... allocate arrays + !------------------------------------------------------------------ + allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr, errmsg=errmsg ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_climo_read_file: tropp_p_in allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate tropp_p_in, error = ' // errmsg) + end if + !------------------------------------------------------------------ + ! ... read in the tropopause pressure + !------------------------------------------------------------------ + ierr = pio_inq_varid( pio_id, 'trop_p', vid ) + start = (/ 1, 1, 1 /) + count = (/ nlon, nlat, ntimes /) + ierr = pio_get_var( pio_id, vid, start, count, tropp_p_in ) + + !------------------------------------------------------------------ + ! ... close the netcdf file + !------------------------------------------------------------------ + call pio_closefile( pio_id ) + + !-------------------------------------------------------------------- + ! ... regrid + !-------------------------------------------------------------------- + + allocate( tropp_p_loc(pcols,ntimes), stat=ierr, errmsg=errmsg ) + + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_climo_read_file: tropp_p_loc allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate tropp_p_loc, error = ' // errmsg) + end if + + call get_rlat_all_p(pcols, to_lats) + call get_rlon_all_p(pcols, to_lons) + call lininterp_init(lon, nlon, to_lons, pcols, 2, lon_wgts, zero, twopi) + call lininterp_init(lat, nlat, to_lats, pcols, 1, lat_wgts) + do n=1,ntimes + call lininterp(tropp_p_in(:,:,n), nlon, nlat, tropp_p_loc(1:pcols,n), pcols, lon_wgts, lat_wgts) + end do + call lininterp_finish(lon_wgts) + call lininterp_finish(lat_wgts) + deallocate(lon) + deallocate(lat) + deallocate(tropp_p_in) + + !-------------------------------------------------------- + ! ... initialize the monthly day of year times + !-------------------------------------------------------- + + allocate( tropp_days(tropp_slices), stat=ierr, errmsg=errmsg ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_climo_read_file: tropp_days allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate tropp_days, error = ' // errmsg) + end if + + do n = 1,tropp_slices + tropp_days(n) = get_calday( dates(n), 0 ) + end do + if (masterproc) then + write(iulog,*) 'tropopause_climo_read_file: tropp_days (fractional day-of-year in climatology dataset) =' + write(iulog,'(1p,5g15.8)') tropp_days(:) + endif + + !-------------------------------------------------------- + ! Mark variables as initialized so they are not read from initial conditions + !-------------------------------------------------------- + call mark_as_initialized('tropopause_air_pressure_from_climatology_dataset') + call mark_as_initialized('tropopause_calendar_days_from_climatology') + call mark_as_initialized('filename_of_tropopause_climatology') + + end subroutine tropopause_climo_read_file +end module tropopause_climo_read diff --git a/src/physics/utils/tropopause_climo_read.meta b/src/physics/utils/tropopause_climo_read.meta new file mode 100644 index 00000000..65dfbcf6 --- /dev/null +++ b/src/physics/utils/tropopause_climo_read.meta @@ -0,0 +1,29 @@ +[ccpp-table-properties] + name = tropopause_climo_read + type = module + +[ccpp-arg-table] + name = tropopause_climo_read + type = module +[ tropp_slices ] + standard_name = number_of_months_in_year + units = 1 + type = integer + dimensions = () +[ tropp_p_loc ] + standard_name = tropopause_air_pressure_from_climatology_dataset + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_months_in_year) +[ tropp_days ] + standard_name = tropopause_calendar_days_from_climatology + long_name = Climatological tropopause calendar day indices from file + units = 1 + type = real | kind = kind_phys + dimensions = (number_of_months_in_year) +[ tropopause_climo_file ] + standard_name = filename_of_tropopause_climatology + long_name = File path to tropopause climatology file + units = none + type = character | kind = len=256 + dimensions = () diff --git a/src/utils/interpolate_data.F90 b/src/utils/interpolate_data.F90 new file mode 100644 index 00000000..c02f5582 --- /dev/null +++ b/src/utils/interpolate_data.F90 @@ -0,0 +1,1230 @@ +module interpolate_data +! Description: +! Routines for interpolation of data in latitude, longitude and time. +! +! Author: Gathered from a number of places and put into the current format by Jim Edwards +! +! Modules Used: +! + use shr_kind_mod, only: r8 => shr_kind_r8 + use cam_abortutils, only: endrun + use cam_logfile, only: iulog + implicit none + private +! +! Public Methods: +! + + public :: interp_type, lininterp, vertinterp, bilin, get_timeinterp_factors + public :: lininterp_init, lininterp_finish + type interp_type + real(r8), pointer :: wgts(:) + real(r8), pointer :: wgtn(:) + integer, pointer :: jjm(:) + integer, pointer :: jjp(:) + end type interp_type + interface lininterp + module procedure lininterp_original + module procedure lininterp_full1d + module procedure lininterp1d + module procedure lininterp2d2d + module procedure lininterp2d1d + module procedure lininterp3d2d + end interface + +integer, parameter, public :: extrap_method_zero = 0 +integer, parameter, public :: extrap_method_bndry = 1 +integer, parameter, public :: extrap_method_cycle = 2 + +contains + subroutine lininterp_full1d (arrin, yin, nin, arrout, yout, nout) + integer, intent(in) :: nin, nout + real(r8), intent(in) :: arrin(nin), yin(nin), yout(nout) + real(r8), intent(out) :: arrout(nout) + type (interp_type) :: interp_wgts + + call lininterp_init(yin, nin, yout, nout, extrap_method_bndry, interp_wgts) + call lininterp1d(arrin, nin, arrout, nout, interp_wgts) + call lininterp_finish(interp_wgts) + + end subroutine lininterp_full1d + + subroutine lininterp_init(yin, nin, yout, nout, extrap_method, interp_wgts, & + cyclicmin, cyclicmax) +! +! Description: +! Initialize a variable of type(interp_type) with weights for linear interpolation. +! this variable can then be used in calls to lininterp1d and lininterp2d. +! yin is a 1d array of length nin of locations to interpolate from - this array must +! be monotonic but can be increasing or decreasing +! yout is a 1d array of length nout of locations to interpolate to, this array need +! not be ordered +! extrap_method determines how to handle yout points beyond the bounds of yin +! if 0 set values outside output grid to 0 +! if 1 set to boundary value +! if 2 set to cyclic boundaries +! optional values cyclicmin and cyclicmax can be used to set the bounds of the +! cyclic mapping - these default to 0 and 360. +! + + integer, intent(in) :: nin + integer, intent(in) :: nout + real(r8), intent(in) :: yin(:) ! input mesh + real(r8), intent(in) :: yout(:) ! output mesh + integer, intent(in) :: extrap_method ! if 0 set values outside output grid to 0 + ! if 1 set to boundary value + ! if 2 set to cyclic boundaries + real(r8), intent(in), optional :: cyclicmin, cyclicmax + + type (interp_type), intent(out) :: interp_wgts + real(r8) :: cmin, cmax + real(r8) :: extrap + real(r8) :: dyinwrap + real(r8) :: ratio + real(r8) :: avgdyin + integer :: i, j, icount + integer :: jj + real(r8), pointer :: wgts(:) + real(r8), pointer :: wgtn(:) + integer, pointer :: jjm(:) + integer, pointer :: jjp(:) + logical :: increasing + ! + ! Check validity of input coordinate arrays: must be monotonically increasing, + ! and have a total of at least 2 elements + ! + if (nin.lt.2) then + call endrun('LININTERP: Must have at least 2 input points for interpolation') + end if + if(present(cyclicmin)) then + cmin=cyclicmin + else + cmin=0_r8 + end if + if(present(cyclicmax)) then + cmax=cyclicmax + else + cmax=360_r8 + end if + if(cmax<=cmin) then + call endrun('LININTERP: cyclic min value must be < max value') + end if + increasing=.true. + icount = 0 + do j=1,nin-1 + if (yin(j).gt.yin(j+1)) icount = icount + 1 + end do + if(icount.eq.nin-1) then + increasing = .false. + icount=0 + endif + if (icount.gt.0) then + call endrun('LININTERP: Non-monotonic input coordinate array found') + end if + allocate(interp_wgts%jjm(nout), & + interp_wgts%jjp(nout), & + interp_wgts%wgts(nout), & + interp_wgts%wgtn(nout)) + + jjm => interp_wgts%jjm + jjp => interp_wgts%jjp + wgts => interp_wgts%wgts + wgtn => interp_wgts%wgtn + + ! + ! Initialize index arrays for later checking + ! + jjm = 0 + jjp = 0 + + extrap = 0._r8 + select case (extrap_method) + case (extrap_method_zero) + ! + ! For values which extend beyond boundaries, set weights + ! such that values will be 0. + ! + do j=1,nout + if(increasing) then + if (yout(j).lt.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).gt.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + else + if (yout(j).gt.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).lt.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + end if + end do + case (extrap_method_bndry) + ! + ! For values which extend beyond boundaries, set weights + ! such that values will just be copied. + ! + do j=1,nout + if(increasing) then + if (yout(j).le.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).gt.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + else + if (yout(j).gt.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).le.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + end if + end do + case (extrap_method_cycle) + ! + ! For values which extend beyond boundaries, set weights + ! for circular boundaries + ! + dyinwrap = yin(1) + (cmax-cmin) - yin(nin) + avgdyin = abs(yin(nin)-yin(1))/(nin-1._r8) + ratio = dyinwrap/avgdyin + if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then + write(iulog,*) 'Lininterp: Bad dyinwrap value =',dyinwrap,& + ' avg=', avgdyin, yin(1),yin(nin) + call endrun('interpolate_data') + end if + + do j=1,nout + if(increasing) then + if (yout(j) <= yin(1)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)-yout(j))/dyinwrap + wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap + else if (yout(j) > yin(nin)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap + wgtn(j) = (yout(j)-yin(nin))/dyinwrap + end if + else + if (yout(j) > yin(1)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)-yout(j))/dyinwrap + wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap + else if (yout(j) <= yin(nin)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap + wgtn(j) = (yout(j)+(cmax-cmin)-yin(nin))/dyinwrap + end if + + endif + end do + end select + + ! + ! Loop though output indices finding input indices and weights + ! + if(increasing) then + do j=1,nout + do jj=1,nin-1 + if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) + wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) + exit + end if + end do + end do + else + do j=1,nout + do jj=1,nin-1 + if (yout(j).le.yin(jj) .and. yout(j).gt.yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) + wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) + exit + end if + end do + end do + end if + + ! + ! Check that interp/extrap points have been found for all outputs + ! + icount = 0 + do j=1,nout + if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1 + ratio=wgts(j)+wgtn(j) + if((ratio<0.9_r8.or.ratio>1.1_r8).and.extrap_method.ne.0) then + write(iulog,*) j, wgts(j),wgtn(j),jjm(j),jjp(j), increasing,extrap_method + call endrun('Bad weight computed in LININTERP_init') + end if + end do + if (icount.gt.0) then + call endrun('LININTERP: Point found without interp indices') + end if + + end subroutine lininterp_init + + subroutine lininterp1d (arrin, n1, arrout, m1, interp_wgts) + !----------------------------------------------------------------------- + ! + ! Purpose: Do a linear interpolation from input mesh to output + ! mesh with weights as set in lininterp_init. + ! + ! + ! Author: Jim Edwards + ! + !----------------------------------------------------------------------- + !----------------------------------------------------------------------- + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1 ! number of input latitudes + integer, intent(in) :: m1 ! number of output latitudes + + real(r8), intent(in) :: arrin(n1) ! input array of values to interpolate + type(interp_type), intent(in) :: interp_wgts + real(r8), intent(out) :: arrout(m1) ! interpolated array + + ! + ! Local workspace + ! + integer j ! latitude indices + integer, pointer :: jjm(:) + integer, pointer :: jjp(:) + + real(r8), pointer :: wgts(:) + real(r8), pointer :: wgtn(:) + + + jjm => interp_wgts%jjm + jjp => interp_wgts%jjp + wgts => interp_wgts%wgts + wgtn => interp_wgts%wgtn + + ! + ! Do the interpolation + ! + do j=1,m1 + arrout(j) = arrin(jjm(j))*wgts(j) + arrin(jjp(j))*wgtn(j) + end do + + return + end subroutine lininterp1d + + subroutine lininterp2d2d(arrin, n1, n2, arrout, m1, m2, wgt1, wgt2) + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1, n2, m1, m2 + real(r8), intent(in) :: arrin(n1,n2) ! input array of values to interpolate + type(interp_type), intent(in) :: wgt1, wgt2 + real(r8), intent(out) :: arrout(m1,m2) ! interpolated array + ! + ! locals + ! + integer i,j ! indices + integer, pointer :: iim(:), jjm(:) + integer, pointer :: iip(:), jjp(:) + + real(r8), pointer :: wgts1(:), wgts2(:) + real(r8), pointer :: wgtn1(:), wgtn2(:) + + real(r8) :: arrtmp(n1,m2) + + + jjm => wgt2%jjm + jjp => wgt2%jjp + wgts2 => wgt2%wgts + wgtn2 => wgt2%wgtn + + iim => wgt1%jjm + iip => wgt1%jjp + wgts1 => wgt1%wgts + wgtn1 => wgt1%wgtn + + do j=1,m2 + do i=1,n1 + arrtmp(i,j) = arrin(i,jjm(j))*wgts2(j) + arrin(i,jjp(j))*wgtn2(j) + end do + end do + + do j=1,m2 + do i=1,m1 + arrout(i,j) = arrtmp(iim(i),j)*wgts1(i) + arrtmp(iip(i),j)*wgtn1(i) + end do + end do + + end subroutine lininterp2d2d + subroutine lininterp2d1d(arrin, n1, n2, arrout, m1, wgt1, wgt2, fldname) + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1, n2, m1 + real(r8), intent(in) :: arrin(n1,n2) ! input array of values to interpolate + type(interp_type), intent(in) :: wgt1, wgt2 + real(r8), intent(out) :: arrout(m1) ! interpolated array + character(len=*), intent(in), optional :: fldname(:) + ! + ! locals + ! + integer i ! indices + integer, pointer :: iim(:), jjm(:) + integer, pointer :: iip(:), jjp(:) + + real(r8), pointer :: wgts(:), wgte(:) + real(r8), pointer :: wgtn(:), wgtw(:) + + jjm => wgt2%jjm + jjp => wgt2%jjp + wgts => wgt2%wgts + wgtn => wgt2%wgtn + + iim => wgt1%jjm + iip => wgt1%jjp + wgtw => wgt1%wgts + wgte => wgt1%wgtn + + do i=1,m1 + arrout(i) = arrin(iim(i),jjm(i))*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i))*wgte(i)*wgts(i) + & + arrin(iim(i),jjp(i))*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i))*wgte(i)*wgtn(i) + end do + + + end subroutine lininterp2d1d + subroutine lininterp3d2d(arrin, n1, n2, n3, arrout, m1, len1, wgt1, wgt2) + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1, n2, n3, m1, len1 ! m1 is to len1 as ncols is to pcols + real(r8), intent(in) :: arrin(n1,n2,n3) ! input array of values to interpolate + type(interp_type), intent(in) :: wgt1, wgt2 + real(r8), intent(out) :: arrout(len1, n3) ! interpolated array + + ! + ! locals + ! + integer i, k ! indices + integer, pointer :: iim(:), jjm(:) + integer, pointer :: iip(:), jjp(:) + + real(r8), pointer :: wgts(:), wgte(:) + real(r8), pointer :: wgtn(:), wgtw(:) + + jjm => wgt2%jjm + jjp => wgt2%jjp + wgts => wgt2%wgts + wgtn => wgt2%wgtn + + iim => wgt1%jjm + iip => wgt1%jjp + wgtw => wgt1%wgts + wgte => wgt1%wgtn + + do k=1,n3 + do i=1,m1 + arrout(i,k) = arrin(iim(i),jjm(i),k)*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i),k)*wgte(i)*wgts(i) + & + arrin(iim(i),jjp(i),k)*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i),k)*wgte(i)*wgtn(i) + end do + end do + + end subroutine lininterp3d2d + + + + + subroutine lininterp_finish(interp_wgts) + type(interp_type) :: interp_wgts + + deallocate(interp_wgts%jjm, & + interp_wgts%jjp, & + interp_wgts%wgts, & + interp_wgts%wgtn) + + nullify(interp_wgts%jjm, & + interp_wgts%jjp, & + interp_wgts%wgts, & + interp_wgts%wgtn) + end subroutine lininterp_finish + + subroutine lininterp_original (arrin, yin, nlev, nlatin, arrout, & + yout, nlatout) + !----------------------------------------------------------------------- + ! + ! Purpose: Do a linear interpolation from input mesh defined by yin to output + ! mesh defined by yout. Where extrapolation is necessary, values will + ! be copied from the extreme edge of the input grid. Vectorization is over + ! the vertical (nlev) dimension. + ! + ! Method: Check validity of input, then determine weights, then do the N-S interpolation. + ! + ! Author: Jim Rosinski + ! Modified: Jim Edwards so that there is no requirement of monotonacity on the yout array + ! + !----------------------------------------------------------------------- + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: nlev ! number of vertical levels + integer, intent(in) :: nlatin ! number of input latitudes + integer, intent(in) :: nlatout ! number of output latitudes + + real(r8), intent(in) :: arrin(nlev,nlatin) ! input array of values to interpolate + real(r8), intent(in) :: yin(nlatin) ! input mesh + real(r8), intent(in) :: yout(nlatout) ! output mesh + + real(r8), intent(out) :: arrout(nlev,nlatout) ! interpolated array + ! + ! Local workspace + ! + integer j, jj ! latitude indices + integer jjprev ! latitude indices + integer k ! level index + integer icount ! number of values + + real(r8) extrap ! percent grid non-overlap + ! + ! Dynamic + ! + integer :: jjm(nlatout) + integer :: jjp(nlatout) + + real(r8) :: wgts(nlatout) + real(r8) :: wgtn(nlatout) + ! + ! Check validity of input coordinate arrays: must be monotonically increasing, + ! and have a total of at least 2 elements + ! + if (nlatin.lt.2) then + call endrun('LININTERP: Must have at least 2 input points for interpolation') + end if + + icount = 0 + do j=1,nlatin-1 + if (yin(j).gt.yin(j+1)) icount = icount + 1 + end do + + + if (icount.gt.0) then + call endrun('LININTERP: Non-monotonic coordinate array(s) found') + end if + ! + ! Initialize index arrays for later checking + ! + do j=1,nlatout + jjm(j) = 0 + jjp(j) = 0 + end do + ! + ! For values which extend beyond N and S boundaries, set weights + ! such that values will just be copied. + ! + extrap = 0._r8 + + do j=1,nlatout + if (yout(j).le.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap=extrap+1._r8 + else if (yout(j).gt.yin(nlatin)) then + jjm(j) = nlatin + jjp(j) = nlatin + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap=extrap+1._r8 + endif + end do + + ! + ! Loop though output indices finding input indices and weights + ! + do j=1,nlatout + do jj=1,nlatin-1 + if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) + wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) + exit + end if + end do + end do + ! + ! Check that interp/extrap points have been found for all outputs + ! + icount = 0 + do j=1,nlatout + if (jjm(j).eq.0 .or. jjp(j).eq.0) then + icount = icount + 1 + end if + end do + if (icount.gt.0) then + call endrun('LININTERP: Point found without interp indices') + end if + ! + ! Do the interpolation + ! + do j=1,nlatout + do k=1,nlev + arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j) + end do + end do + + return + end subroutine lininterp_original + + + subroutine bilin (arrin, xin, yin, nlondin, nlonin, & + nlevdin, nlev, nlatin, arrout, xout, & + yout, nlondout, nlonout, nlevdout, nlatout) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! + ! Do a bilinear interpolation from input mesh defined by xin, yin to output + ! mesh defined by xout, yout. Circularity is assumed in the x-direction so + ! input x-grid must be in degrees east and must start from Greenwich. When + ! extrapolation is necessary in the N-S direction, values will be copied + ! from the extreme edge of the input grid. Vectorization is over the + ! longitude dimension. Input grid is assumed rectangular. Output grid + ! is assumed ragged, i.e. xout is a 2-d mesh. + ! + ! Method: Interpolate first in longitude, then in latitude. + ! + ! Author: Jim Rosinski + ! + !----------------------------------------------------------------------- + use shr_kind_mod, only: r8 => shr_kind_r8 + use cam_abortutils, only: endrun + !----------------------------------------------------------------------- + implicit none + !----------------------------------------------------------------------- + ! + ! Input arguments + ! + integer, intent(in) :: nlondin ! longitude dimension of input grid + integer, intent(in) :: nlonin ! number of real longitudes (input) + integer, intent(in) :: nlevdin ! vertical dimension of input grid + integer, intent(in) :: nlev ! number of vertical levels + integer, intent(in) :: nlatin ! number of input latitudes + integer, intent(in) :: nlatout ! number of output latitudes + integer, intent(in) :: nlondout ! longitude dimension of output grid + integer, intent(in) :: nlonout(nlatout) ! number of output longitudes per lat + integer, intent(in) :: nlevdout ! vertical dimension of output grid + + real(r8), intent(in) :: arrin(nlondin,nlevdin,nlatin) ! input array of values to interpolate + real(r8), intent(in) :: xin(nlondin) ! input x mesh + real(r8), intent(in) :: yin(nlatin) ! input y mesh + real(r8), intent(in) :: xout(nlondout,nlatout) ! output x mesh + real(r8), intent(in) :: yout(nlatout) ! output y mesh + ! + ! Output arguments + ! + real(r8), intent(out) :: arrout(nlondout,nlevdout,nlatout) ! interpolated array + ! + ! Local workspace + ! + integer :: i, ii, iw, ie, iiprev ! longitude indices + integer :: j, jj, js, jn, jjprev ! latitude indices + integer :: k ! level index + integer :: icount ! number of bad values + + real(r8) :: extrap ! percent grid non-overlap + real(r8) :: dxinwrap ! delta-x on input grid for 2-pi + real(r8) :: avgdxin ! avg input delta-x + real(r8) :: ratio ! compare dxinwrap to avgdxin + real(r8) :: sum ! sum of weights (used for testing) + ! + ! Dynamic + ! + integer :: iim(nlondout) ! interpolation index to the left + integer :: iip(nlondout) ! interpolation index to the right + integer :: jjm(nlatout) ! interpolation index to the south + integer :: jjp(nlatout) ! interpolation index to the north + + real(r8) :: wgts(nlatout) ! interpolation weight to the north + real(r8) :: wgtn(nlatout) ! interpolation weight to the north + real(r8) :: wgte(nlondout) ! interpolation weight to the north + real(r8) :: wgtw(nlondout) ! interpolation weight to the north + real(r8) :: igrid(nlonin) ! interpolation weight to the north + ! + ! Check validity of input coordinate arrays: must be monotonically increasing, + ! and have a total of at least 2 elements + ! + if (nlonin < 2 .or. nlatin < 2) then + call endrun ('BILIN: Must have at least 2 input points for interpolation') + end if + + if (xin(1) < 0._r8 .or. xin(nlonin) > 360._r8) then + call endrun ('BILIN: Input x-grid must be between 0 and 360') + end if + + icount = 0 + do i=1,nlonin-1 + if (xin(i) >= xin(i+1)) icount = icount + 1 + end do + + do j=1,nlatin-1 + if (yin(j) >= yin(j+1)) icount = icount + 1 + end do + + do j=1,nlatout-1 + if (yout(j) >= yout(j+1)) icount = icount + 1 + end do + + do j=1,nlatout + do i=1,nlonout(j)-1 + if (xout(i,j) >= xout(i+1,j)) icount = icount + 1 + end do + end do + + if (icount > 0) then + call endrun ('BILIN: Non-monotonic coordinate array(s) found') + end if + + if (yout(nlatout) <= yin(1) .or. yout(1) >= yin(nlatin)) then + call endrun ('BILIN: No overlap in y-coordinates') + end if + + do j=1,nlatout + if (xout(1,j) < 0._r8 .or. xout(nlonout(j),j) > 360._r8) then + call endrun ('BILIN: Output x-grid must be between 0 and 360') + end if + + if (xout(nlonout(j),j) <= xin(1) .or. & + xout(1,j) >= xin(nlonin)) then + call endrun ('BILIN: No overlap in x-coordinates') + end if + end do + ! + ! Initialize index arrays for later checking + ! + do j=1,nlatout + jjm(j) = 0 + jjp(j) = 0 + end do + ! + ! For values which extend beyond N and S boundaries, set weights + ! such that values will just be copied. + ! + do js=1,nlatout + if (yout(js) > yin(1)) exit + jjm(js) = 1 + jjp(js) = 1 + wgts(js) = 1._r8 + wgtn(js) = 0._r8 + end do + + do jn=nlatout,1,-1 + if (yout(jn) <= yin(nlatin)) exit + jjm(jn) = nlatin + jjp(jn) = nlatin + wgts(jn) = 1._r8 + wgtn(jn) = 0._r8 + end do + ! + ! Loop though output indices finding input indices and weights + ! + jjprev = 1 + do j=js,jn + do jj=jjprev,nlatin-1 + if (yout(j) > yin(jj) .and. yout(j) <= yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1) - yout(j)) / (yin(jj+1) - yin(jj)) + wgtn(j) = (yout(j) - yin(jj)) / (yin(jj+1) - yin(jj)) + goto 30 + end if + end do + call endrun ('BILIN: Failed to find interp values') +30 jjprev = jj + end do + + dxinwrap = xin(1) + 360._r8 - xin(nlonin) + ! + ! Check for sane dxinwrap values. Allow to differ no more than 10% from avg + ! + avgdxin = (xin(nlonin)-xin(1))/(nlonin-1._r8) + ratio = dxinwrap/avgdxin + if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then + write(iulog,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin + call endrun('BILIN: Insane dxinwrap value') + end if + ! + ! Check that interp/extrap points have been found for all outputs, and that + ! they are reasonable (i.e. within 32-bit roundoff) + ! + icount = 0 + do j=1,nlatout + if (jjm(j) == 0 .or. jjp(j) == 0) icount = icount + 1 + sum = wgts(j) + wgtn(j) + if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1 + if (wgts(j) < 0._r8 .or. wgts(j) > 1._r8) icount = icount + 1 + if (wgtn(j) < 0._r8 .or. wgtn(j) > 1._r8) icount = icount + 1 + end do + + if (icount > 0) then + call endrun ('BILIN: something bad in latitude indices or weights') + end if + ! + ! Do the bilinear interpolation + ! + do j=1,nlatout + ! + ! Initialize index arrays for later checking + ! + do i=1,nlondout + iim(i) = 0 + iip(i) = 0 + end do + ! + ! For values which extend beyond E and W boundaries, set weights such that + ! values will be interpolated between E and W edges of input grid. + ! + do iw=1,nlonout(j) + if (xout(iw,j) > xin(1)) exit + iim(iw) = nlonin + iip(iw) = 1 + wgtw(iw) = (xin(1) - xout(iw,j)) /dxinwrap + wgte(iw) = (xout(iw,j)+360._r8 - xin(nlonin))/dxinwrap + end do + + do ie=nlonout(j),1,-1 + if (xout(ie,j) <= xin(nlonin)) exit + iim(ie) = nlonin + iip(ie) = 1 + wgtw(ie) = (xin(1)+360._r8 - xout(ie,j)) /dxinwrap + wgte(ie) = (xout(ie,j) - xin(nlonin))/dxinwrap + end do + ! + ! Loop though output indices finding input indices and weights + ! + iiprev = 1 + do i=iw,ie + do ii=iiprev,nlonin-1 + if (xout(i,j) > xin(ii) .and. xout(i,j) <= xin(ii+1)) then + iim(i) = ii + iip(i) = ii + 1 + wgtw(i) = (xin(ii+1) - xout(i,j)) / (xin(ii+1) - xin(ii)) + wgte(i) = (xout(i,j) - xin(ii)) / (xin(ii+1) - xin(ii)) + goto 60 + end if + end do + call endrun ('BILIN: Failed to find interp values') +60 iiprev = ii + end do + + icount = 0 + do i=1,nlonout(j) + if (iim(i) == 0 .or. iip(i) == 0) icount = icount + 1 + sum = wgtw(i) + wgte(i) + if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1 + if (wgtw(i) < 0._r8 .or. wgtw(i) > 1._r8) icount = icount + 1 + if (wgte(i) < 0._r8 .or. wgte(i) > 1._r8) icount = icount + 1 + end do + + if (icount > 0) then + write(iulog,*)'BILIN: j=',j,' Something bad in longitude indices or weights' + call endrun('BILIN: Something bad in longitude indices or weights') + end if + ! + ! Do the interpolation, 1st in longitude then latitude + ! + do k=1,nlev + do i=1,nlonin + igrid(i) = arrin(i,k,jjm(j))*wgts(j) + arrin(i,k,jjp(j))*wgtn(j) + end do + + do i=1,nlonout(j) + arrout(i,k,j) = igrid(iim(i))*wgtw(i) + igrid(iip(i))*wgte(i) + end do + end do + end do + + + return + end subroutine bilin + +!========================================================================================= + +subroutine vertinterp(ncol, ncold, nlev, pin, pout, arrin, arrout, & + extrapolate, ln_interp, ps, phis, tbot) + + ! Vertically interpolate input array to output pressure level. The + ! interpolation is linear in either pressure (the default) or ln pressure. + ! + ! If above the top boundary then return top boundary value. + ! + ! If below bottom boundary then use bottom boundary value, or optionally + ! for T or Z use the extrapolation algorithm from ECMWF (which was taken + ! from the NCL code base). + + use physconst, only: rair, rga + + !------------------------------Arguments-------------------------------- + integer, intent(in) :: ncol ! number active columns + integer, intent(in) :: ncold ! column dimension + integer, intent(in) :: nlev ! vertical dimension + real(r8), intent(in) :: pin(ncold,nlev) ! input pressure levels + real(r8), intent(in) :: pout ! output pressure level + real(r8), intent(in) :: arrin(ncold,nlev) ! input array + real(r8), intent(out) :: arrout(ncold) ! output array (interpolated) + + character(len=1), optional, intent(in) :: extrapolate ! either 'T' or 'Z' + logical, optional, intent(in) :: ln_interp ! set true to interolate + ! in ln(pressure) + real(r8), optional, intent(in) :: ps(ncold) ! surface pressure + real(r8), optional, intent(in) :: phis(ncold) ! surface geopotential + real(r8), optional, intent(in) :: tbot(ncold) ! temperature at bottom level + + !---------------------------Local variables----------------------------- + real(r8) :: alpha + logical :: linear_interp + logical :: do_extrapolate_T, do_extrapolate_Z + + integer :: i,k ! indices + integer :: kupper(ncold) ! Level indices for interpolation + real(r8) :: dpu ! upper level pressure difference + real(r8) :: dpl ! lower level pressure difference + logical :: found(ncold) ! true if input levels found + logical :: error ! error flag + !---------------------------------------------------------------------------- + + alpha = 0.0065_r8*rair*rga + + do_extrapolate_T = .false. + do_extrapolate_Z = .false. + if (present(extrapolate)) then + + if (extrapolate == 'T') do_extrapolate_T = .true. + if (extrapolate == 'Z') do_extrapolate_Z = .true. + + if (.not. do_extrapolate_T .and. .not. do_extrapolate_Z) then + call endrun ('VERTINTERP: extrapolate must be T or Z') + end if + if (.not. present(ps)) then + call endrun ('VERTINTERP: ps required for extrapolation') + end if + if (.not. present(phis)) then + call endrun ('VERTINTERP: phis required for extrapolation') + end if + if (do_extrapolate_Z) then + if (.not. present(tbot)) then + call endrun ('VERTINTERP: extrapolate must be T or Z') + end if + end if + end if + + linear_interp = .true. + if (present(ln_interp)) then + if (ln_interp) linear_interp = .false. + end if + + do i = 1, ncol + found(i) = .false. + kupper(i) = 1 + end do + error = .false. + + ! Find indices of upper pressure bound + do k = 1, nlev - 1 + do i = 1, ncol + if ((.not. found(i)) .and. pin(i,k)= pin(i,nlev)) then + + if (do_extrapolate_T) then + ! use ECMWF algorithm to extrapolate T + arrout(i) = extrapolate_T() + + else if (do_extrapolate_Z) then + ! use ECMWF algorithm to extrapolate Z + arrout(i) = extrapolate_Z() + + else + ! use bottom value + arrout(i) = arrin(i,nlev) + end if + + else if (found(i)) then + if (linear_interp) then + ! linear interpolation in p + dpu = pout - pin(i,kupper(i)) + dpl = pin(i,kupper(i)+1) - pout + arrout(i) = (arrin(i,kupper(i))*dpl + arrin(i,kupper(i)+1)*dpu)/(dpl + dpu) + else + ! linear interpolation in ln(p) + arrout(i) = arrin(i,kupper(i)) + (arrin(i,kupper(i)+1) - arrin(i,kupper(i)))* & + log(pout/pin(i,kupper(i))) / & + log(pin(i,kupper(i)+1)/pin(i,kupper(i))) + end if + else + error = .true. + end if + end do + + ! Error check + if (error) then + call endrun ('VERTINTERP: ERROR FLAG') + end if + +contains + + !====================================================================================== + + real(r8) function extrapolate_T() + + ! local variables + real(r8) :: sixth + real(r8) :: tstar + real(r8) :: hgt + real(r8) :: alnp + real(r8) :: t0 + real(r8) :: tplat + real(r8) :: tprime0 + !------------------------------------------------------------------------- + + sixth = 1._r8 / 6._r8 + tstar = arrin(i,nlev)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8)) + hgt = phis(i)*rga + + if (hgt < 2000._r8) then + alnp = alpha*log(pout/ps(i)) + else + t0 = tstar + 0.0065_r8*hgt + tplat = min(t0, 298._r8) + + if (hgt <= 2500._r8) then + tprime0 = 0.002_r8*((2500._r8 - hgt)*t0 + & + (hgt - 2000._r8)*tplat) + else + tprime0 = tplat + end if + + if (tprime0 < tstar) then + alnp = 0._r8 + else + alnp = rair*(tprime0 - tstar)/phis(i) * log(pout/ps(i)) + end if + + end if + + extrapolate_T = tstar*(1._r8 + alnp + 0.5_r8*alnp**2 + sixth*alnp**3) + + end function extrapolate_T + + !====================================================================================== + + real(r8) function extrapolate_Z() + + ! local variables + real(r8) :: sixth + real(r8) :: hgt + real(r8) :: tstar + real(r8) :: t0 + real(r8) :: alph + real(r8) :: alnp + !------------------------------------------------------------------------- + + sixth = 1._r8 / 6._r8 + hgt = phis(i)*rga + tstar = tbot(i)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8)) + t0 = tstar + 0.0065_r8*hgt + + if (tstar <= 290.5_r8 .and. t0 > 290.5_r8) then + alph = rair/phis(i) * (290.5_r8 - tstar) + + else if (tstar > 290.5_r8 .and. t0 > 290.5_r8) then + alph = 0._r8 + tstar = 0.5_r8*(290.5_r8 + tstar) + + else + alph = alpha + + end if + + if (tstar < 255._r8) then + tstar = 0.5_r8*(tstar + 255._r8) + end if + alnp = alph*log(pout/ps(i)) + + extrapolate_Z = hgt - rair*tstar*rga*log(pout/ps(i))* & + (1._r8 + 0.5_r8*alnp + sixth*alnp**2) + + end function extrapolate_Z + +end subroutine vertinterp + +!========================================================================================= + + subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, & + fact1, fact2, str) + !--------------------------------------------------------------------------- + ! + ! Purpose: Determine time interpolation factors (normally for a boundary dataset) + ! for linear interpolation. + ! + ! Method: Assume 365 days per year. Output variable fact1 will be the weight to + ! apply to data at calendar time "cdayminus", and fact2 the weight to apply + ! to data at time "cdayplus". Combining these values will produce a result + ! valid at time "cday". Output arguments fact1 and fact2 will be between + ! 0 and 1, and fact1 + fact2 = 1 to roundoff. + ! + ! Author: Jim Rosinski + ! + !--------------------------------------------------------------------------- + implicit none + ! + ! Arguments + ! + logical, intent(in) :: cycflag ! flag indicates whether dataset is being cycled yearly + + integer, intent(in) :: np1 ! index points to forward time slice matching cdayplus + + real(r8), intent(in) :: cdayminus ! calendar day of rearward time slice + real(r8), intent(in) :: cdayplus ! calendar day of forward time slice + real(r8), intent(in) :: cday ! calenar day to be interpolated to + real(r8), intent(out) :: fact1 ! time interpolation factor to apply to rearward time slice + real(r8), intent(out) :: fact2 ! time interpolation factor to apply to forward time slice + + character(len=*), intent(in) :: str ! string to be added to print in case of error (normally the callers name) + ! + ! Local workspace + ! + real(r8) :: deltat ! time difference (days) between cdayminus and cdayplus + real(r8), parameter :: daysperyear = 365._r8 ! number of days in a year + ! + ! Initial sanity checks + ! + if (np1 == 1 .and. .not. cycflag) then + call endrun ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed') + end if + + if (np1 < 1) then + call endrun ('GETFACTORS:'//str//' input arg np1 must be > 0') + end if + + if (cycflag) then + if ((cday < 1._r8) .or. (cday > (daysperyear+1._r8))) then + write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday + call endrun ('get_timeinterp_factors GETFACTORS bad cday (1)') + end if + else + if (cday < 1._r8) then + write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday + call endrun ('get_timeinterp_factors GETFACTORS bad cday (2)') + end if + end if + ! + ! Determine time interpolation factors. Account for December-January + ! interpolation if dataset is being cycled yearly. + ! + if (cycflag .and. np1 == 1) then ! Dec-Jan interpolation + deltat = cdayplus + daysperyear - cdayminus + if (cday > cdayplus) then ! We are in December + fact1 = (cdayplus + daysperyear - cday)/deltat + fact2 = (cday - cdayminus)/deltat + else ! We are in January + fact1 = (cdayplus - cday)/deltat + fact2 = (cday + daysperyear - cdayminus)/deltat + end if + else + deltat = cdayplus - cdayminus + fact1 = (cdayplus - cday)/deltat + fact2 = (cday - cdayminus)/deltat + end if + + if (.not. valid_timeinterp_factors (fact1, fact2)) then + write(iulog,*) 'GETFACTORS: ', str, ' bad fact1 and/or fact2=', fact1, fact2 + call endrun ('get_timeinterp_factors GETFACTORS bad fact1 and/or fact2') + end if + + return + end subroutine get_timeinterp_factors + + logical function valid_timeinterp_factors (fact1, fact2) + !--------------------------------------------------------------------------- + ! + ! Purpose: check sanity of time interpolation factors to within 32-bit roundoff + ! + !--------------------------------------------------------------------------- + implicit none + + real(r8), intent(in) :: fact1, fact2 ! time interpolation factors + + valid_timeinterp_factors = .true. + + ! The fact1 .ne. fact1 and fact2 .ne. fact2 comparisons are to detect NaNs. + if (abs(fact1+fact2-1._r8) > 1.e-6_r8 .or. & + fact1 > 1.000001_r8 .or. fact1 < -1.e-6_r8 .or. & + fact2 > 1.000001_r8 .or. fact2 < -1.e-6_r8 .or. & + fact1 .ne. fact1 .or. fact2 .ne. fact2) then + + valid_timeinterp_factors = .false. + end if + + return + end function valid_timeinterp_factors + +end module interpolate_data