mom_horizontal_regridding Module Reference

Detailed Description

Horizontal interpolation.

Data Types
interface	fill_boundaries
	Fill grid edges. More...
interface	horiz_interp_and_extrap_tracer
	Extrapolate and interpolate data. More...

Functions/Subroutines
subroutine	fill_miss_2d (aout, good, fill, prev, G, smooth, num_pass, relc, crit, keep_bug, debug)
	Use ICE-9 algorithm to populate points (fill=1) with valid data (good=1). If no information is available, Then use a previous guess (prev). Optionally (smooth) blend the filled points to achieve a more desirable result. More...
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, public	mystats (array, missing, is, ie, js, je, k, mesg)
	Write to the terminal some basic statistics about the k-th level of an array. 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, spongeOngrid, m_to_Z)
	Extrapolate and interpolate using a FMS time interpolation handle. More...
subroutine	meshgrid (x, y, x_T, y_T)
	Create a 2d-mesh of grid coordinates from 1-d arrays. More...
integer function, dimension(0:size(m, 1)+1, 0:size(m, 2)+1)	fill_boundaries_int (m, cyclic_x, tripolar_n)
	Fill grid edges for integer data. More...
real function, dimension(0:size(m, 1)+1, 0:size(m, 2)+1)	fill_boundaries_real (m, cyclic_x, tripolar_n)
	Fill grid edges for real data. More...
subroutine	smooth_heights (zi, fill, bad, sor, niter, cyclic_x, tripolar_n)
	Solve del2 (zi) = 0 using successive iterations with a 5 point stencil. Only points fill==1 are modified. Except where bad==1, information propagates isotropically in index space. The resulting solution in each region is an approximation to del2(zi)=0 subject to boundary conditions along the valid points curve bounding this region. More...

Function/Subroutine Documentation

fill_boundaries_int()

integer function, dimension(0:size(m,1)+1,0:size(m,2)+1) mom_horizontal_regridding::fill_boundaries_int ( integer, dimension(:,:), intent(in)  m, logical, intent(in)  cyclic_x, logical, intent(in)  tripolar_n  )

private

Fill grid edges for integer data.

Parameters

[in]	m	input array (ND)
[in]	cyclic_x	True if domain is zonally re-entrant
[in]	tripolar_n	True if domain has an Arctic fold

Definition at line 900 of file MOM_horizontal_regridding.F90.

  integer, dimension(:,:), intent(in)             :: m !< input array (ND)
  logical,                 intent(in)             :: cyclic_x !< True if domain is zonally re-entrant
  logical,                 intent(in)             :: tripolar_n !< True if domain has an Arctic fold
  integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp
 
  real,    dimension(size(m,1),size(m,2))         :: m_real
  real,    dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real
 
  m_real = real(m)
 
  mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n)
 
  mp = int(mp_real)
 

References fill_boundaries_real().

Here is the call graph for this function:

fill_boundaries_real()

real function, dimension(0:size(m,1)+1,0:size(m,2)+1) mom_horizontal_regridding::fill_boundaries_real ( real, dimension(:,:), intent(in)  m, logical, intent(in)  cyclic_x, logical, intent(in)  tripolar_n  )

private

Fill grid edges for real data.

Parameters

[in]	m	input array (ND)
[in]	cyclic_x	True if domain is zonally re-entrant
[in]	tripolar_n	True if domain has an Arctic fold

Definition at line 918 of file MOM_horizontal_regridding.F90.

  real, dimension(:,:), intent(in)             :: m !< input array (ND)
  logical,              intent(in)             :: cyclic_x !< True if domain is zonally re-entrant
  logical,              intent(in)             :: tripolar_n !< True if domain has an Arctic fold
  real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp
 
  integer :: ni,nj,i,j
 
  ni=size(m,1); nj=size(m,2)
 
  mp(1:ni,1:nj)=m(:,:)
 
  if (cyclic_x) then
    mp(0,1:nj)=m(ni,1:nj)
    mp(ni+1,1:nj)=m(1,1:nj)
  else
    mp(0,1:nj)=m(1,1:nj)
    mp(ni+1,1:nj)=m(ni,1:nj)
  endif
 
  mp(1:ni,0)=m(1:ni,1)
  if (tripolar_n) then
    do i=1,ni
      mp(i,nj+1)=m(ni-i+1,nj)
    enddo
  else
    mp(1:ni,nj+1)=m(1:ni,nj)
  endif
 

Referenced by fill_boundaries_int().

Here is the caller graph for this function:

fill_miss_2d()

subroutine mom_horizontal_regridding::fill_miss_2d ( real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)  aout, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  good, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  fill, real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional  prev, type(ocean_grid_type), intent(inout)  G, logical, intent(in), optional  smooth, integer, intent(in), optional  num_pass, real, intent(in), optional  relc, real, intent(in), optional  crit, logical, intent(in), optional  keep_bug, logical, intent(in), optional  debug  )

private

Use ICE-9 algorithm to populate points (fill=1) with valid data (good=1). If no information is available, Then use a previous guess (prev). Optionally (smooth) blend the filled points to achieve a more desirable result.

Parameters

[in,out]	g	The ocean's grid structure.
[in,out]	aout	The array with missing values to fill
[in]	good	Valid data mask for incoming array
[in]	fill	Same shape array of points which need
[in]	prev	First guess where isolated holes exist.
[in]	smooth	If present and true, apply a number of Laplacian iterations to the interpolated data
[in]	num_pass	The maximum number of iterations
[in]	relc	A relaxation coefficient for Laplacian (ND)
[in]	crit	A minimal value for deltas between iterations.
[in]	keep_bug	Use an algorithm with a bug that dates to the "sienna" code release.
[in]	debug	If true, write verbose debugging messages.

Definition at line 104 of file MOM_horizontal_regridding.F90.

  use mom_coms, only : sum_across_pes
 
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)), &
                         intent(inout) :: aout !< The array with missing values to fill
  real, dimension(SZI_(G),SZJ_(G)), &
                         intent(in)    :: good !< Valid data mask for incoming array
                                               !! (1==good data; 0==missing data).
  real, dimension(SZI_(G),SZJ_(G)), &
                         intent(in)    :: fill !< Same shape array of points which need
                                               !! filling (1==fill;0==dont fill)
  real, dimension(SZI_(G),SZJ_(G)), &
               optional, intent(in)    :: prev !< First guess where isolated holes exist.
  logical,     optional, intent(in)    :: smooth !< If present and true, apply a number of
                                                 !! Laplacian iterations to the interpolated data
  integer,     optional, intent(in)    :: num_pass !< The maximum number of iterations
  real,        optional, intent(in)    :: relc !< A relaxation coefficient for Laplacian (ND)
  real,        optional, intent(in)    :: crit !< A minimal value for deltas between iterations.
  logical,     optional, intent(in)    :: keep_bug !< Use an algorithm with a bug that dates
                                                   !! to the "sienna" code release.
  logical,     optional, intent(in)    :: debug !< If true, write verbose debugging messages.
 
  real, dimension(SZI_(G),SZJ_(G)) :: b,r
  real, dimension(SZI_(G),SZJ_(G)) :: fill_pts, good_, good_new
 
  character(len=256) :: mesg  ! The text of an error message
  integer :: i,j,k
  real    :: east,west,north,south,sor
  real    :: ge,gw,gn,gs,ngood
  logical :: do_smooth,siena_bug
  real    :: nfill, nfill_prev
  integer, parameter :: num_pass_default = 10000
  real, parameter :: relc_default = 0.25, crit_default = 1.e-3
 
  integer :: npass
  integer :: is, ie, js, je
  real    :: relax_coeff, acrit, ares
  logical :: debug_it
 
  debug_it=.false.
  if (PRESENT(debug)) debug_it=debug
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
 
  npass = num_pass_default
  if (PRESENT(num_pass)) npass = num_pass
 
  relax_coeff = relc_default
  if (PRESENT(relc)) relax_coeff = relc
 
  acrit = crit_default
  if (PRESENT(crit)) acrit = crit
 
  siena_bug=.false.
  if (PRESENT(keep_bug)) siena_bug = keep_bug
 
  do_smooth=.false.
  if (PRESENT(smooth)) do_smooth=smooth
 
  fill_pts(:,:) = fill(:,:)
 
  nfill = sum(fill(is:ie,js:je))
  call sum_across_pes(nfill)
 
  nfill_prev = nfill
  good_(:,:) = good(:,:)
  r(:,:) = 0.0
 
  do while (nfill > 0.0)
 
    call pass_var(good_,g%Domain)
    call pass_var(aout,g%Domain)
 
    b(:,:)=aout(:,:)
    good_new(:,:)=good_(:,:)
 
    do j=js,je ; do i=is,ie
 
      if (good_(i,j) == 1.0 .or. fill(i,j) == 0.) cycle
 
      ge=good_(i+1,j) ; gw=good_(i-1,j)
      gn=good_(i,j+1) ; gs=good_(i,j-1)
      east=0.0 ; west=0.0 ; north=0.0 ; south=0.0
      if (ge == 1.0) east = aout(i+1,j)*ge
      if (gw == 1.0) west = aout(i-1,j)*gw
      if (gn == 1.0) north = aout(i,j+1)*gn
      if (gs == 1.0) south = aout(i,j-1)*gs
 
      ngood = ge+gw+gn+gs
      if (ngood > 0.) then
        b(i,j)=(east+west+north+south)/ngood
        !### Replace this with
        ! b(i,j) = ((east+west) + (north+south))/ngood
        fill_pts(i,j) = 0.0
        good_new(i,j) = 1.0
      endif
    enddo ; enddo
 
    aout(is:ie,js:je) = b(is:ie,js:je)
    good_(is:ie,js:je) = good_new(is:ie,js:je)
    nfill_prev = nfill
    nfill = sum(fill_pts(is:ie,js:je))
    call sum_across_pes(nfill)
 
    if (nfill == nfill_prev .and. PRESENT(prev)) then
      do j=js,je ; do i=is,ie ; if (fill_pts(i,j) == 1.0) then
        aout(i,j) = prev(i,j)
        fill_pts(i,j) = 0.0
      endif ; enddo ; enddo
    elseif (nfill == nfill_prev) then
      call mom_error(warning, &
           'Unable to fill missing points using either data at the same vertical level from a connected basin'//&
           'or using a point from a previous vertical level.  Make sure that the original data has some valid'//&
           'data in all basins.', .true.)
      write(mesg,*) 'nfill=',nfill
      call mom_error(warning, mesg, .true.)
    endif
 
    nfill = sum(fill_pts(is:ie,js:je))
    call sum_across_pes(nfill)
 
  enddo
 
  if (do_smooth) then ; do k=1,npass
    call pass_var(aout,g%Domain)
    do j=js,je ; do i=is,ie
      if (fill(i,j) == 1) then
        east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j))
        north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1))
        r(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + &
                              west*aout(i-1,j)+east*aout(i+1,j) - &
                             (south+north+west+east)*aout(i,j))
        !### Appropriate parentheses should be added here, but they will change answers.
        ! r(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + &
        !                         (west*aout(i-1,j)+east*aout(i+1,j))) - &
        !                        ((south+north)+(west+east))*aout(i,j) )
      else
        r(i,j) = 0.
      endif
    enddo ; enddo
    ares = 0.0
    do j=js,je ; do i=is,ie
      aout(i,j) = r(i,j) + aout(i,j)
      ares = max(ares, abs(r(i,j)))
    enddo ; enddo
    call max_across_pes(ares)
    if (ares <= acrit) exit
  enddo ; endif
 
  do j=js,je ; do i=is,ie
    if (good_(i,j) == 0.0 .and. fill_pts(i,j) == 1.0) then
      write(mesg,*) 'In fill_miss, fill, good,i,j= ',fill_pts(i,j),good_(i,j),i,j
      call mom_error(warning, mesg, .true.)
      call mom_error(fatal,"MOM_initialize: "// &
           "fill is true and good is false after fill_miss, how did this happen? ")
    endif
 enddo ; enddo
 

Referenced by horiz_interp_and_extrap_tracer_fms_id(), and horiz_interp_and_extrap_tracer_record().

Here is the caller graph for this function:

horiz_interp_and_extrap_tracer_fms_id()

subroutine mom_horizontal_regridding::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, logical, intent(in), optional  spongeOngrid, 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]	spongeongrid	If present and true, the sponge data are on the model grid
[in]	m_to_z	A conversion factor from meters to the units of depth. If missing, GbathyT must be in m.

Definition at line 591 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
  logical,     optional, intent(in)    :: spongeOngrid !< If present and true, the sponge data are on the model grid
  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.
  logical :: spongeDataOngrid
  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)
 
  spongedataongrid=.false.
  if (PRESENT(spongeongrid)) spongedataongrid=spongeongrid
  if (.not. spongedataongrid) then
    allocate(lon_in(id),lat_in(jd))
    call mpp_get_axis_data(axes_data(1), lon_in)
    call mpp_get_axis_data(axes_data(2), lat_in)
  endif
 
  allocate(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(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)
 
  if (.not. spongedataongrid) then
    ! 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
    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
  else
    allocate(data_in(isd:ied,jsd:jed,kd))
  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)
 
 
  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
 
  !  roundoff = 3.0*EPSILON(missing_value)
  roundoff = 1.e-4
 
  if (.not.spongedataongrid) then
    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
    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=.true.)
      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
  else
      call time_interp_external(fms_id, time, data_in, verbose=.true.)
      do k=1,kd
        do j=js,je
          do i=is,ie
            tr_z(i,j,k)=data_in(i,j,k)
            if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
          enddo
        enddo
      enddo
  endif
 

References fill_miss_2d(), meshgrid(), and mystats().

Here is the call graph for this function:

horiz_interp_and_extrap_tracer_record()

subroutine mom_horizontal_regridding::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
	)

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 268 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
 
  roundoff = 3.0*epsilon(missing_value)
 
  ! loop through each data level and interpolate to model grid.
  ! after interpolating, fill in points which will be needed
  ! to define the layers
  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=.true.)
    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
 

References fill_miss_2d(), meshgrid(), and mystats().

Here is the call graph for this function:

meshgrid()

subroutine mom_horizontal_regridding::meshgrid ( real, dimension(:), intent(in)  x, real, dimension(:), intent(in)  y, real, dimension(size(x,1),size(y,1)), intent(inout)  x_T, real, dimension(size(x,1),size(y,1)), intent(inout)  y_T  )

private

Create a 2d-mesh of grid coordinates from 1-d arrays.

Parameters

[in]	x	input 1-dimensional vector
[in]	y	input 1-dimensional vector
[in,out]	x_t	output 2-dimensional array
[in,out]	y_t	output 2-dimensional array

Definition at line 877 of file MOM_horizontal_regridding.F90.

  real, dimension(:),                   intent(in)    :: x  !< input 1-dimensional vector
  real, dimension(:),                   intent(in)    :: y  !< input 1-dimensional vector
  real, dimension(size(x,1),size(y,1)), intent(inout) :: x_T !< output 2-dimensional array
  real, dimension(size(x,1),size(y,1)), intent(inout) :: y_T !< output 2-dimensional array
 
  integer :: ni,nj,i,j
 
  ni=size(x,1) ; nj=size(y,1)
 
  do j=1,nj ; do i=1,ni
    x_t(i,j) = x(i)
  enddo ; enddo
 
  do j=1,nj ; do i=1,ni
    y_t(i,j) = y(j)
  enddo ; enddo
 

Referenced by horiz_interp_and_extrap_tracer_fms_id(), and horiz_interp_and_extrap_tracer_record().

Here is the caller graph for this function:

mystats()

subroutine, public mom_horizontal_regridding::mystats	(	real, dimension(:,:), intent(in)	array,
		real, intent(in)	missing,
		integer	is,
		integer	ie,
		integer	js,
		integer	je,
		integer	k,
		character(len=*)	mesg
	)

Write to the terminal some basic statistics about the k-th level of an array.

Parameters

[in]	array	input array (ND)
[in]	missing	missing value (ND)
	is	Horizontal loop bounds to calculate statistics for
	ie	Horizontal loop bounds to calculate statistics for
	js	Horizontal loop bounds to calculate statistics for
	je	Horizontal loop bounds to calculate statistics for
	k	Level to calculate statistics for
	mesg	Label to use in message

Definition at line 61 of file MOM_horizontal_regridding.F90.

  real, dimension(:,:), intent(in) :: array !< input array (ND)
  real, intent(in) :: missing !< missing value (ND)
  !!@{
  !> Horizontal loop bounds to calculate statistics for
  integer :: is,ie,js,je
  !!@}
  integer :: k !< Level to calculate statistics for
  character(len=*) :: mesg !< Label to use in message
  ! Local variables
  real :: minA, maxA
  integer :: i,j
  logical :: found
  character(len=120) :: lMesg
  mina = 9.e24 ; maxa = -9.e24 ; found = .false.
 
  do j=js,je ; do i=is,ie
    if (array(i,j) /= array(i,j)) stop 'Nan!'
    if (abs(array(i,j)-missing) > 1.e-6*abs(missing)) then
      if (found) then
        mina = min(mina, array(i,j))
        maxa = max(maxa, array(i,j))
      else
        found = .true.
        mina = array(i,j)
        maxa = array(i,j)
      endif
    endif
  enddo ; enddo
  call min_across_pes(mina)
  call max_across_pes(maxa)
  if (is_root_pe()) then
    write(lmesg(1:120),'(2(a,es12.4),a,i3,x,a)') &
         'init_from_Z: min=',mina,' max=',maxa,' Level=',k,trim(mesg)
    call mom_mesg(lmesg,2)
  endif
 

Referenced by horiz_interp_and_extrap_tracer_fms_id(), horiz_interp_and_extrap_tracer_record(), and mom_tracer_initialization_from_z::mom_initialize_tracer_from_z().

Here is the caller graph for this function:

smooth_heights()

subroutine mom_horizontal_regridding::smooth_heights ( real, dimension(:,:), intent(inout)  zi, integer, dimension(size(zi,1),size(zi,2)), intent(in)  fill, integer, dimension(size(zi,1),size(zi,2)), intent(in)  bad, real, intent(in)  sor, integer, intent(in)  niter, logical, intent(in)  cyclic_x, logical, intent(in)  tripolar_n  )

private

Solve del2 (zi) = 0 using successive iterations with a 5 point stencil. Only points fill==1 are modified. Except where bad==1, information propagates isotropically in index space. The resulting solution in each region is an approximation to del2(zi)=0 subject to boundary conditions along the valid points curve bounding this region.

Parameters

[in,out]	zi	input and output array (ND)
[in]	fill	same shape as zi, 1=fill
[in]	bad	same shape as zi, 1=bad data
[in]	sor	relaxation coefficient (ND)
[in]	niter	maximum number of iterations
[in]	cyclic_x	true if domain is zonally reentrant
[in]	tripolar_n	true if domain has an Arctic fold

Definition at line 955 of file MOM_horizontal_regridding.F90.

  real,    dimension(:,:),                   intent(inout) :: zi !< input and output array (ND)
  integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< same shape as zi, 1=fill
  integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad  !< same shape as zi, 1=bad data
  real,                                      intent(in)  :: sor !< relaxation coefficient (ND)
  integer,                                   intent(in) :: niter !< maximum number of iterations
  logical,                                   intent(in) :: cyclic_x !< true if domain is zonally reentrant
  logical,                                   intent(in) :: tripolar_n !< true if domain has an Arctic fold
 
  ! Local variables
  real, dimension(size(zi,1),size(zi,2)) :: res, m
  integer, dimension(size(zi,1),size(zi,2),4) :: B
  real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp
  integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm
  integer :: i,j,k,n
  integer :: ni,nj
  real :: Isum, bsum
 
  ni=size(zi,1) ; nj=size(zi,2)
 
 
  mp(:,:) = fill_boundaries(zi,cyclic_x,tripolar_n)
 
  b(:,:,:) = 0.0
  nm(:,:) = fill_boundaries(bad,cyclic_x,tripolar_n)
 
  do j=1,nj
    do i=1,ni
      if (fill(i,j) == 1) then
        b(i,j,1)=1-nm(i+1,j);b(i,j,2)=1-nm(i-1,j)
        b(i,j,3)=1-nm(i,j+1);b(i,j,4)=1-nm(i,j-1)
      endif
    enddo
  enddo
 
  do n=1,niter
    do j=1,nj
      do i=1,ni
        if (fill(i,j) == 1) then
          bsum = real(b(i,j,1)+b(i,j,2)+b(i,j,3)+b(i,j,4))
          isum = 1.0/bsum
          res(i,j)=isum*(b(i,j,1)*mp(i+1,j)+b(i,j,2)*mp(i-1,j)+&
                   b(i,j,3)*mp(i,j+1)+b(i,j,4)*mp(i,j-1)) - mp(i,j)
        endif
      enddo
    enddo
    res(:,:)=res(:,:)*sor
 
    do j=1,nj
      do i=1,ni
        mp(i,j)=mp(i,j)+res(i,j)
      enddo
    enddo
 
    zi(:,:)=mp(1:ni,1:nj)
    mp = fill_boundaries(zi,cyclic_x,tripolar_n)
  enddo