Detailed Description

Extrapolate and interpolate data.

Definition at line 52 of file MOM_horizontal_regridding.F90.

Private functions
subroutine	horiz_interp_and_extrap_tracer_record (filename, varnam, conversion, recnum, G, tr_z, mask_z, z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, homogenize, m_to_Z)
	Extrapolate and interpolate from a file record. More...

subroutine	horiz_interp_and_extrap_tracer_fms_id (fms_id, Time, conversion, G, tr_z, mask_z, z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, homogenize, m_to_Z)
	Extrapolate and interpolate using a FMS time interpolation handle. More...

Functions and subroutines

◆ horiz_interp_and_extrap_tracer_fms_id()

subroutine mom_horizontal_regridding::horiz_interp_and_extrap_tracer::horiz_interp_and_extrap_tracer_fms_id	(	integer, intent(in)	fms_id,
		type(time_type), intent(in)	Time,
		real, intent(in)	conversion,
		type(ocean_grid_type), intent(inout)	G,
		real, dimension(:,:,:), allocatable	tr_z,
		real, dimension(:,:,:), allocatable	mask_z,
		real, dimension(:), allocatable	z_in,
		real, dimension(:), allocatable	z_edges_in,
		real, intent(out)	missing_value,
		logical, intent(in)	reentrant_x,
		logical, intent(in)	tripolar_n,
		logical, intent(in), optional	homogenize,
		real, intent(in), optional	m_to_Z
	)

private

Extrapolate and interpolate using a FMS time interpolation handle.

Parameters

[in]	fms_id	A unique id used by the FMS time interpolator
[in]	time	A FMS time type
[in]	conversion	Conversion factor for tracer.
[in,out]	g	Grid object
	tr_z	pointer to allocatable tracer array on local model grid and native vertical levels.
	mask_z	pointer to allocatable tracer mask array on local model grid and native vertical levels.
	z_in	Cell grid values for input data.
	z_edges_in	Cell grid edge values for input data. (Intent out)
[out]	missing_value	The missing value in the returned array.
[in]	reentrant_x	If true, this grid is reentrant in the x-direction
[in]	tripolar_n	If true, this is a northern tripolar grid
[in]	homogenize	If present and true, horizontally homogenize data to produce perfectly "flat" initial conditions
[in]	m_to_z	A conversion factor from meters to the units of depth. If missing, GbathyT must be in m.

Definition at line 604 of file MOM_horizontal_regridding.F90.

  
   integer,               intent(in)    :: fms_id     !< A unique id used by the FMS time interpolator
   type(time_type),       intent(in)    :: Time       !< A FMS time type
   real,                  intent(in)    :: conversion !< Conversion factor for tracer.
   type(ocean_grid_type), intent(inout) :: G          !< Grid object
   real, allocatable, dimension(:,:,:)  :: tr_z       !< pointer to allocatable tracer array on local
                                                      !! model grid and native vertical levels.
   real, allocatable, dimension(:,:,:)  :: mask_z     !< pointer to allocatable tracer mask array on
                                                      !! local model grid and native vertical levels.
   real, allocatable,     dimension(:)  :: z_in       !< Cell grid values for input data.
   real, allocatable,     dimension(:)  :: z_edges_in !< Cell grid edge values for input data. (Intent out)
   real,                  intent(out)   :: missing_value !< The missing value in the returned array.
   logical,               intent(in)    :: reentrant_x !< If true, this grid is reentrant in the x-direction
   logical,               intent(in)    :: tripolar_n !< If true, this is a northern tripolar grid
   logical,     optional, intent(in)    :: homogenize !< If present and true, horizontally homogenize data
                                                      !! to produce perfectly "flat" initial conditions
   real,        optional, intent(in)    :: m_to_Z     !< A conversion factor from meters to the units
                                                      !! of depth.  If missing, G%bathyT must be in m.
  
   ! Local variables
   real, dimension(:,:),  allocatable   :: tr_in,tr_inp !< A 2-d array for holding input data on
                                                      !! native horizontal grid and extended grid
                                                      !! with poles.
   real, dimension(:,:,:), allocatable  :: data_in    !< A buffer for storing the full 3-d time-interpolated array
                                                      !! on the original grid
   real, dimension(:,:),  allocatable   :: mask_in    !< A 2-d mask for extended input grid.
  
   real :: PI_180
   integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp
   integer :: i,j,k
   integer, dimension(4) :: start, count, dims, dim_id
   real, dimension(:,:), allocatable :: x_in, y_in
   real, dimension(:), allocatable  :: lon_in, lat_in
   real, dimension(:), allocatable  :: lat_inp, last_row
   real :: max_lat, min_lat, pole, max_depth, npole
   real :: roundoff  ! The magnitude of roundoff, usually ~2e-16.
   logical :: add_np
   character(len=8)  :: laynum
   type(horiz_interp_type) :: Interp
   type(axistype), dimension(4) :: axes_data
   integer :: is, ie, js, je     ! compute domain indices
   integer :: isc,iec,jsc,jec    ! global compute domain indices
   integer :: isg, ieg, jsg, jeg ! global extent
   integer :: isd, ied, jsd, jed ! data domain indices
   integer :: id_clock_read
   integer, dimension(4) :: fld_sz
   character(len=12)  :: dim_name(4)
   logical :: debug=.false.
   real :: npoints,varAvg
   real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out
   real, dimension(SZI_(G),SZJ_(G)) :: good, fill
   real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev
   real, dimension(SZI_(G),SZJ_(G))  :: good2,fill2
   real, dimension(SZI_(G),SZJ_(G))  :: nlevs
  
   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
   isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
   isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
  
   id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
  
  
   pi_180=atan(1.0)/45.
  
   ! Open NetCDF file and if present, extract data and spatial coordinate information
   ! The convention adopted here requires that the data be written in (i,j,k) ordering.
  
   call cpu_clock_begin(id_clock_read)
  
   fld_sz = get_external_field_size(fms_id)
  
   if (allocated(lon_in)) deallocate(lon_in)
   if (allocated(lat_in)) deallocate(lat_in)
   if (allocated(z_in)) deallocate(z_in)
   if (allocated(z_edges_in)) deallocate(z_edges_in)
   if (allocated(tr_z)) deallocate(tr_z)
   if (allocated(mask_z)) deallocate(mask_z)
  
   axes_data =  get_external_field_axes(fms_id)
  
   id = fld_sz(1) ; jd  = fld_sz(2) ; kd = fld_sz(3)
   allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1))
   allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))
  
   call mpp_get_axis_data(axes_data(1), lon_in)
   call mpp_get_axis_data(axes_data(2), lat_in)
   call mpp_get_axis_data(axes_data(3), z_in)
  
   if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
  
   call cpu_clock_end(id_clock_read)
  
   missing_value = get_external_field_missing(fms_id)
  
  
 ! extrapolate the input data to the north pole using the northerm-most latitude
  
   max_lat = maxval(lat_in)
   add_np=.false.
   if (max_lat < 90.0) then
     add_np = .true.
     jdp = jd+1
     allocate(lat_inp(jdp))
     lat_inp(1:jd) = lat_in(:)
     lat_inp(jd+1) = 90.0
     deallocate(lat_in)
     allocate(lat_in(1:jdp))
     lat_in(:) = lat_inp(:)
   else
     jdp=jd
   endif
  
 ! construct level cell boundaries as the mid-point between adjacent centers
  
   z_edges_in(1) = 0.0
   do k=2,kd
    z_edges_in(k) = 0.5*(z_in(k-1)+z_in(k))
   enddo
   z_edges_in(kd+1) = 2.0*z_in(kd) - z_in(kd-1)
  
   call horiz_interp_init()
  
   lon_in = lon_in*pi_180
   lat_in = lat_in*pi_180
   allocate(x_in(id,jdp), y_in(id,jdp))
   call meshgrid(lon_in, lat_in, x_in, y_in)
  
   lon_out(:,:) = g%geoLonT(:,:)*pi_180
   lat_out(:,:) = g%geoLatT(:,:)*pi_180
  
   allocate(data_in(id,jd,kd)) ; data_in(:,:,:)=0.0
   allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0
   allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0
   allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0
   allocate(last_row(id))    ; last_row(:)=0.0
  
   max_depth = maxval(g%bathyT)
   call mpp_max(max_depth)
  
   if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
  
   if (is_root_pe()) &
   call time_interp_external(fms_id, time, data_in, verbose=.true.)
  
 ! loop through each data level and interpolate to model grid.
 ! after interpolating, fill in points which will be needed
 ! to define the layers
  
 !  roundoff = 3.0*EPSILON(missing_value)
   roundoff = 1.e-4
  
   do k=1,kd
     write(laynum,'(I8)') k ; laynum = adjustl(laynum)
  
     if (is_root_pe()) then
       tr_in(1:id,1:jd) = data_in(1:id,1:jd,k)
       if (add_np) then
          last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0
          do i=1,id
             if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then
                pole = pole+last_row(i)
                npole = npole+1.0
             endif
          enddo
          if (npole > 0) then
             pole=pole/npole
          else
             pole=missing_value
          endif
          tr_inp(:,1:jd) = tr_in(:,:)
          tr_inp(:,jdp) = pole
       else
          tr_inp(:,:) = tr_in(:,:)
       endif
  
     endif
  
     call mpp_sync()
     call mpp_broadcast(tr_inp, id*jdp, root_pe())
     call mpp_sync_self()
  
     mask_in=0.0
  
     do j=1,jdp ; do i=1,id
       if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then
         mask_in(i,j)=1.0
         tr_inp(i,j) = tr_inp(i,j) * conversion
       else
         tr_inp(i,j) = missing_value
       endif
     enddo ; enddo
  
     ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
     if (k == 1) then
       call horiz_interp_new(interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
                interp_method='bilinear', src_modulo=reentrant_x)
     endif
  
     if (debug) then
       call mystats(tr_in, missing_value, 1, id, 1, jd, k, 'Tracer from file')
     endif
  
     tr_out(:,:) = 0.0
  
     call horiz_interp(interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, &
                       new_missing_handle=.true.)
  
     mask_out(:,:) = 1.0
     do j=js,je ; do i=is,ie
       if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j) = 0.
     enddo ; enddo
  
     fill(:,:) = 0.0 ; good(:,:) = 0.0
  
     npoints = 0 ; varavg = 0.
     do j=js,je ; do i=is,ie
       if (mask_out(i,j) < 1.0) then
         tr_out(i,j) = missing_value
       else
         good(i,j) = 1.0
         npoints = npoints + 1
         varavg = varavg + tr_out(i,j)
       endif
       if ((g%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= g%bathyT(i,j)) .and. &
           (mask_out(i,j) < 1.0)) &
         fill(i,j)=1.0
     enddo ;  enddo
     call pass_var(fill, g%Domain)
     call pass_var(good, g%Domain)
  
     if (debug) then
       call mystats(tr_out, missing_value, is, ie, js, je, k, 'variable from horiz_interp()')
     endif
  
     ! Horizontally homogenize data to produce perfectly "flat" initial conditions
     if (PRESENT(homogenize)) then ; if (homogenize) then
       call sum_across_pes(npoints)
       call sum_across_pes(varavg)
       if (npoints>0) then
          varavg = varavg/real(npoints)
       endif
       tr_out(:,:) = varavg
     endif ; endif
  
   ! tr_out contains input z-space data on the model grid with missing values
   ! now fill in missing values using "ICE-nine" algorithm.
  
     tr_outf(:,:) = tr_out(:,:)
     if (k==1) tr_prev(:,:) = tr_outf(:,:)
     good2(:,:) = good(:,:)
     fill2(:,:) = fill(:,:)
  
     call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, smooth=.true.)
  
 !    if (debug) then
 !      call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI)
 !    endif
  
 !    call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()')
  
     tr_z(:,:,k) = tr_outf(:,:)*g%mask2dT(:,:)
     mask_z(:,:,k) = good2(:,:) + fill2(:,:)
     tr_prev(:,:) = tr_z(:,:,k)
  
     if (debug) then
       call hchksum(tr_prev,'field after fill ',g%HI)
     endif
  
   enddo ! kd
  

◆ horiz_interp_and_extrap_tracer_record()

subroutine mom_horizontal_regridding::horiz_interp_and_extrap_tracer::horiz_interp_and_extrap_tracer_record	(	character(len=*), intent(in)	filename,
		character(len=*), intent(in)	varnam,
		real, intent(in)	conversion,
		integer, intent(in)	recnum,
		type(ocean_grid_type), intent(inout)	G,
		real, dimension(:,:,:), allocatable	tr_z,
		real, dimension(:,:,:), allocatable	mask_z,
		real, dimension(:), allocatable	z_in,
		real, dimension(:), allocatable	z_edges_in,
		real, intent(out)	missing_value,
		logical, intent(in)	reentrant_x,
		logical, intent(in)	tripolar_n,
		logical, intent(in), optional	homogenize,
		real, intent(in), optional	m_to_Z
	)

private

Extrapolate and interpolate from a file record.

Parameters

[in]	filename	Path to file containing tracer to be interpolated.
[in]	varnam	Name of tracer in filee.
[in]	conversion	Conversion factor for tracer.
[in]	recnum	Record number of tracer to be read.
[in,out]	g	Grid object
	tr_z	pointer to allocatable tracer array on local model grid and input-file vertical levels.
	mask_z	pointer to allocatable tracer mask array on local model grid and input-file vertical levels.
	z_in	Cell grid values for input data.
	z_edges_in	Cell grid edge values for input data.
[out]	missing_value	The missing value in the returned array.
[in]	reentrant_x	If true, this grid is reentrant in the x-direction
[in]	tripolar_n	If true, this is a northern tripolar grid
[in]	homogenize	If present and true, horizontally homogenize data to produce perfectly "flat" initial conditions
[in]	m_to_z	A conversion factor from meters to the units of depth. If missing, GbathyT must be in m.

Definition at line 270 of file MOM_horizontal_regridding.F90.

  
   character(len=*),      intent(in)    :: filename   !< Path to file containing tracer to be
                                                      !! interpolated.
   character(len=*),      intent(in)    :: varnam     !< Name of tracer in filee.
   real,                  intent(in)    :: conversion !< Conversion factor for tracer.
   integer,               intent(in)    :: recnum     !< Record number of tracer to be read.
   type(ocean_grid_type), intent(inout) :: G          !< Grid object
   real, allocatable, dimension(:,:,:)  :: tr_z       !< pointer to allocatable tracer array on local
                                                      !! model grid and input-file vertical levels.
   real, allocatable, dimension(:,:,:)  :: mask_z     !< pointer to allocatable tracer mask array on
                                                      !! local model grid and input-file vertical levels.
   real, allocatable,     dimension(:)  :: z_in       !< Cell grid values for input data.
   real, allocatable,     dimension(:)  :: z_edges_in !< Cell grid edge values for input data.
   real,                  intent(out)   :: missing_value !< The missing value in the returned array.
   logical,               intent(in)    :: reentrant_x !< If true, this grid is reentrant in the x-direction
   logical,               intent(in)    :: tripolar_n !< If true, this is a northern tripolar grid
   logical,     optional, intent(in)    :: homogenize !< If present and true, horizontally homogenize data
                                                      !! to produce perfectly "flat" initial conditions
   real,        optional, intent(in)    :: m_to_Z     !< A conversion factor from meters to the units
                                                      !! of depth.  If missing, G%bathyT must be in m.
  
   ! Local variables
   real, dimension(:,:),  allocatable   :: tr_in, tr_inp ! A 2-d array for holding input data on
                                                      ! native horizontal grid and extended grid
                                                      ! with poles.
   real, dimension(:,:),  allocatable   :: mask_in    ! A 2-d mask for extended input grid.
  
   real :: PI_180
   integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp
   integer :: i,j,k
   integer, dimension(4) :: start, count, dims, dim_id
   real, dimension(:,:), allocatable :: x_in, y_in
   real, dimension(:), allocatable  :: lon_in, lat_in
   real, dimension(:), allocatable  :: lat_inp, last_row
   real :: max_lat, min_lat, pole, max_depth, npole
   real :: roundoff  ! The magnitude of roundoff, usually ~2e-16.
   real :: add_offset, scale_factor
   logical :: add_np
   character(len=8)  :: laynum
   type(horiz_interp_type) :: Interp
   integer :: is, ie, js, je     ! compute domain indices
   integer :: isc,iec,jsc,jec    ! global compute domain indices
   integer :: isg, ieg, jsg, jeg ! global extent
   integer :: isd, ied, jsd, jed ! data domain indices
   integer :: id_clock_read
   character(len=12)  :: dim_name(4)
   logical :: debug=.false.
   real :: npoints,varAvg
   real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out
   real, dimension(SZI_(G),SZJ_(G)) :: good, fill
   real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev
   real, dimension(SZI_(G),SZJ_(G))  :: good2,fill2
   real, dimension(SZI_(G),SZJ_(G))  :: nlevs
  
   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
   isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
   isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
  
   id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
  
  
   if (allocated(tr_z)) deallocate(tr_z)
   if (allocated(mask_z)) deallocate(mask_z)
   if (allocated(z_edges_in)) deallocate(z_edges_in)
  
   pi_180=atan(1.0)/45.
  
   ! Open NetCDF file and if present, extract data and spatial coordinate information
   ! The convention adopted here requires that the data be written in (i,j,k) ordering.
  
   call cpu_clock_begin(id_clock_read)
  
  
   rcode = nf90_open(filename, nf90_nowrite, ncid)
   if (rcode /= 0) call mom_error(fatal,"error opening file "//trim(filename)//&
                            " in hinterp_extrap")
   rcode = nf90_inq_varid(ncid, varnam, varid)
   if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(varnam)//&
                                  " in file "//trim(filename)//" in hinterp_extrap")
  
   rcode = nf90_inquire_variable(ncid, varid, ndims=ndims, dimids=dims)
   if (rcode /= 0) call mom_error(fatal,'error inquiring dimensions hinterp_extrap')
   if (ndims < 3) call mom_error(fatal,"Variable "//trim(varnam)//" in file "// &
               trim(filename)//" has too few dimensions.")
  
   rcode = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)
   if (rcode /= 0) call mom_error(fatal,"error reading dimension 1 data for "// &
                 trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
   rcode = nf90_inq_varid(ncid, dim_name(1), dim_id(1))
   if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(1))//&
                                  " in file "//trim(filename)//" in hinterp_extrap")
   rcode = nf90_inquire_dimension(ncid, dims(2), dim_name(2), len=jd)
   if (rcode /= 0) call mom_error(fatal,"error reading dimension 2 data for "// &
                 trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
   rcode = nf90_inq_varid(ncid, dim_name(2), dim_id(2))
   if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(2))//&
                                  " in file "//trim(filename)//" in hinterp_extrap")
   rcode = nf90_inquire_dimension(ncid, dims(3), dim_name(3), len=kd)
   if (rcode /= 0) call mom_error(fatal,"error reading dimension 3 data for "// &
                 trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
   rcode = nf90_inq_varid(ncid, dim_name(3), dim_id(3))
   if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(3))//&
                                  " in file "//trim(filename)//" in hinterp_extrap")
  
  
   missing_value=0.0
   rcode = nf90_get_att(ncid, varid, "_FillValue", missing_value)
   if (rcode /= 0) call mom_error(fatal,"error finding missing value for "//&
        trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
  
   rcode = nf90_get_att(ncid, varid, "add_offset", add_offset)
   if (rcode /= 0) add_offset = 0.0
  
   rcode = nf90_get_att(ncid, varid, "scale_factor", scale_factor)
   if (rcode /= 0) scale_factor = 1.0
  
  
   if (allocated(lon_in)) deallocate(lon_in)
   if (allocated(lat_in)) deallocate(lat_in)
   if (allocated(z_in)) deallocate(z_in)
   if (allocated(z_edges_in)) deallocate(z_edges_in)
   if (allocated(tr_z)) deallocate(tr_z)
   if (allocated(mask_z)) deallocate(mask_z)
  
  
   allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1))
   allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))
  
   start = 1; count = 1; count(1) = id
   rcode = nf90_get_var(ncid, dim_id(1), lon_in, start, count)
   if (rcode /= 0) call mom_error(fatal,"error reading dimension 1 values for var_name "// &
                 trim(varnam)//",dim_name "//trim(dim_name(1))//" in file "// trim(filename)//" in hinterp_extrap")
   start = 1; count = 1; count(1) = jd
   rcode = nf90_get_var(ncid, dim_id(2), lat_in, start, count)
   if (rcode /= 0) call mom_error(fatal,"error reading dimension 2 values for var_name "// &
                 trim(varnam)//",dim_name "//trim(dim_name(2))//" in file "// trim(filename)//" in  hinterp_extrap")
   start = 1; count = 1; count(1) = kd
   rcode = nf90_get_var(ncid, dim_id(3), z_in, start, count)
   if (rcode /= 0) call mom_error(fatal,"error reading dimension 3 values for var_name "// &
                 trim(varnam//",dim_name "//trim(dim_name(3)))//" in file "// trim(filename)//" in  hinterp_extrap")
  
   call cpu_clock_end(id_clock_read)
  
   if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
  
 ! extrapolate the input data to the north pole using the northerm-most latitude
  
   max_lat = maxval(lat_in)
   add_np=.false.
   if (max_lat < 90.0) then
     add_np=.true.
     jdp=jd+1
     allocate(lat_inp(jdp))
     lat_inp(1:jd)=lat_in(:)
     lat_inp(jd+1)=90.0
     deallocate(lat_in)
     allocate(lat_in(1:jdp))
     lat_in(:)=lat_inp(:)
   else
     jdp=jd
   endif
  
 ! construct level cell boundaries as the mid-point between adjacent centers
  
   z_edges_in(1) = 0.0
   do k=2,kd
    z_edges_in(k)=0.5*(z_in(k-1)+z_in(k))
   enddo
   z_edges_in(kd+1)=2.0*z_in(kd) - z_in(kd-1)
  
   call horiz_interp_init()
  
   lon_in = lon_in*pi_180
   lat_in = lat_in*pi_180
   allocate(x_in(id,jdp),y_in(id,jdp))
   call meshgrid(lon_in,lat_in, x_in, y_in)
  
   lon_out(:,:) = g%geoLonT(:,:)*pi_180
   lat_out(:,:) = g%geoLatT(:,:)*pi_180
  
  
   allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0
   allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0
   allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0
   allocate(last_row(id))    ; last_row(:)=0.0
  
   max_depth = maxval(g%bathyT)
   call mpp_max(max_depth)
  
   if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
  
  
 ! loop through each data level and interpolate to model grid.
 ! after interpolating, fill in points which will be needed
 ! to define the layers
  
   roundoff = 3.0*epsilon(missing_value)
  
  
   do k=1,kd
     write(laynum,'(I8)') k ; laynum = adjustl(laynum)
  
     if (is_root_pe()) then
       start = 1 ; start(3) = k ; count(:) = 1 ; count(1) = id ; count(2) = jd
       rcode = nf90_get_var(ncid,varid, tr_in, start, count)
       if (rcode /= 0) call mom_error(fatal,"hinterp_and_extract_from_Fie: "//&
            "error reading level "//trim(laynum)//" of variable "//&
            trim(varnam)//" in file "// trim(filename))
  
       if (add_np) then
          last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0
          do i=1,id
             if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then
                pole = pole+last_row(i)
                npole = npole+1.0
             endif
          enddo
          if (npole > 0) then
             pole=pole/npole
          else
             pole=missing_value
          endif
          tr_inp(:,1:jd) = tr_in(:,:)
          tr_inp(:,jdp) = pole
       else
          tr_inp(:,:) = tr_in(:,:)
       endif
  
     endif
  
     call mpp_sync()
     call mpp_broadcast(tr_inp, id*jdp, root_pe())
     call mpp_sync_self()
  
     mask_in=0.0
  
     do j=1,jdp
       do i=1,id
          if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then
            mask_in(i,j) = 1.0
            tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion
          else
            tr_inp(i,j) = missing_value
          endif
       enddo
     enddo
  
  
 ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
  
  
     if (k == 1) then
       call horiz_interp_new(interp,x_in,y_in,lon_out(is:ie,js:je),lat_out(is:ie,js:je), &
                interp_method='bilinear',src_modulo=reentrant_x)
     endif
  
     if (debug) then
        call mystats(tr_inp,missing_value, is,ie,js,je,k,'Tracer from file')
     endif
  
     tr_out(:,:) = 0.0
  
     call horiz_interp(interp,tr_inp,tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.)
  
     mask_out=1.0
     do j=js,je
       do i=is,ie
         if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j)=0.
       enddo
     enddo
  
     fill = 0.0; good = 0.0
  
     npoints = 0 ; varavg = 0.
     do j=js,je
       do i=is,ie
         if (mask_out(i,j) < 1.0) then
           tr_out(i,j)=missing_value
         else
           good(i,j)=1.0
           npoints = npoints + 1
           varavg = varavg + tr_out(i,j)
         endif
         if (g%mask2dT(i,j) == 1.0 .and. z_edges_in(k) <= g%bathyT(i,j) .and. mask_out(i,j) < 1.0) &
           fill(i,j)=1.0
       enddo
     enddo
     call pass_var(fill,g%Domain)
     call pass_var(good,g%Domain)
  
     if (debug) then
       call mystats(tr_out,missing_value, is,ie,js,je,k,'variable from horiz_interp()')
     endif
  
     ! Horizontally homogenize data to produce perfectly "flat" initial conditions
     if (PRESENT(homogenize)) then
        if (homogenize) then
           call sum_across_pes(npoints)
           call sum_across_pes(varavg)
           if (npoints>0) then
              varavg = varavg/real(npoints)
           endif
           tr_out(:,:) = varavg
        endif
     endif
  
 ! tr_out contains input z-space data on the model grid with missing values
 ! now fill in missing values using "ICE-nine" algorithm.
  
     tr_outf(:,:) = tr_out(:,:)
     if (k==1) tr_prev(:,:) = tr_outf(:,:)
     good2(:,:) = good(:,:)
     fill2(:,:) = fill(:,:)
  
     call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, smooth=.true.)
     call mystats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()')
  
     tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
     mask_z(:,:,k) = good2(:,:) + fill2(:,:)
  
     tr_prev(:,:) = tr_z(:,:,k)
  
     if (debug) then
       call hchksum(tr_prev,'field after fill ',g%HI)
     endif
  
   enddo ! kd
  

The documentation for this interface was generated from the following file:

/home/docs/checkouts/readthedocs.org/user_builds/mom6-doctesting/checkouts/brokenbreathe/src/framework/MOM_horizontal_regridding.F90

Detailed Description

Private functions

Functions and subroutines

◆ horiz_interp_and_extrap_tracer_fms_id()

◆ horiz_interp_and_extrap_tracer_record()