MOM6
Private functions | List of all members
mom_horizontal_regridding::horiz_interp_and_extrap_tracer Interface Reference

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, spongeOngrid, 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,
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_idA unique id used by the FMS time interpolator
[in]timeA FMS time type
[in]conversionConversion factor for tracer.
[in,out]gGrid object
tr_zpointer to allocatable tracer array on local model grid and native vertical levels.
mask_zpointer to allocatable tracer mask array on local model grid and native vertical levels.
z_inCell grid values for input data.
z_edges_inCell grid edge values for input data. (Intent out)
[out]missing_valueThe missing value in the returned array.
[in]reentrant_xIf true, this grid is reentrant in the x-direction
[in]tripolar_nIf true, this is a northern tripolar grid
[in]homogenizeIf present and true, horizontally homogenize data to produce perfectly "flat" initial conditions
[in]spongeongridIf present and true, the sponge data are on the model grid
[in]m_to_zA 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.

591 
592  integer, intent(in) :: fms_id !< A unique id used by the FMS time interpolator
593  type(time_type), intent(in) :: Time !< A FMS time type
594  real, intent(in) :: conversion !< Conversion factor for tracer.
595  type(ocean_grid_type), intent(inout) :: G !< Grid object
596  real, allocatable, dimension(:,:,:) :: tr_z !< pointer to allocatable tracer array on local
597  !! model grid and native vertical levels.
598  real, allocatable, dimension(:,:,:) :: mask_z !< pointer to allocatable tracer mask array on
599  !! local model grid and native vertical levels.
600  real, allocatable, dimension(:) :: z_in !< Cell grid values for input data.
601  real, allocatable, dimension(:) :: z_edges_in !< Cell grid edge values for input data. (Intent out)
602  real, intent(out) :: missing_value !< The missing value in the returned array.
603  logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction
604  logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid
605  logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
606  !! to produce perfectly "flat" initial conditions
607  logical, optional, intent(in) :: spongeOngrid !< If present and true, the sponge data are on the model grid
608  real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units
609  !! of depth. If missing, G%bathyT must be in m.
610 
611  ! Local variables
612  real, dimension(:,:), allocatable :: tr_in,tr_inp !< A 2-d array for holding input data on
613  !! native horizontal grid and extended grid
614  !! with poles.
615  real, dimension(:,:,:), allocatable :: data_in !< A buffer for storing the full 3-d time-interpolated array
616  !! on the original grid
617  real, dimension(:,:), allocatable :: mask_in !< A 2-d mask for extended input grid.
618 
619  real :: PI_180
620  integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp
621  integer :: i,j,k
622  integer, dimension(4) :: start, count, dims, dim_id
623  real, dimension(:,:), allocatable :: x_in, y_in
624  real, dimension(:), allocatable :: lon_in, lat_in
625  real, dimension(:), allocatable :: lat_inp, last_row
626  real :: max_lat, min_lat, pole, max_depth, npole
627  real :: roundoff ! The magnitude of roundoff, usually ~2e-16.
628  logical :: add_np
629  character(len=8) :: laynum
630  type(horiz_interp_type) :: Interp
631  type(axistype), dimension(4) :: axes_data
632  integer :: is, ie, js, je ! compute domain indices
633  integer :: isc,iec,jsc,jec ! global compute domain indices
634  integer :: isg, ieg, jsg, jeg ! global extent
635  integer :: isd, ied, jsd, jed ! data domain indices
636  integer :: id_clock_read
637  integer, dimension(4) :: fld_sz
638  character(len=12) :: dim_name(4)
639  logical :: debug=.false.
640  logical :: spongeDataOngrid
641  real :: npoints,varAvg
642  real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out
643  real, dimension(SZI_(G),SZJ_(G)) :: good, fill
644  real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev
645  real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2
646  real, dimension(SZI_(G),SZJ_(G)) :: nlevs
647 
648  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
649  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
650  isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
651 
652  id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
653 
654  pi_180=atan(1.0)/45.
655 
656  ! Open NetCDF file and if present, extract data and spatial coordinate information
657  ! The convention adopted here requires that the data be written in (i,j,k) ordering.
658 
659  call cpu_clock_begin(id_clock_read)
660 
661  fld_sz = get_external_field_size(fms_id)
662  if (allocated(lon_in)) deallocate(lon_in)
663  if (allocated(lat_in)) deallocate(lat_in)
664  if (allocated(z_in)) deallocate(z_in)
665  if (allocated(z_edges_in)) deallocate(z_edges_in)
666  if (allocated(tr_z)) deallocate(tr_z)
667  if (allocated(mask_z)) deallocate(mask_z)
668 
669  axes_data = get_external_field_axes(fms_id)
670 
671  id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3)
672 
673  spongedataongrid=.false.
674  if (PRESENT(spongeongrid)) spongedataongrid=spongeongrid
675  if (.not. spongedataongrid) then
676  allocate(lon_in(id),lat_in(jd))
677  call mpp_get_axis_data(axes_data(1), lon_in)
678  call mpp_get_axis_data(axes_data(2), lat_in)
679  endif
680 
681  allocate(z_in(kd),z_edges_in(kd+1))
682 
683  allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))
684 
685  call mpp_get_axis_data(axes_data(3), z_in)
686 
687  if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
688 
689  call cpu_clock_end(id_clock_read)
690 
691  missing_value = get_external_field_missing(fms_id)
692 
693  if (.not. spongedataongrid) then
694  ! extrapolate the input data to the north pole using the northerm-most latitude
695  max_lat = maxval(lat_in)
696  add_np=.false.
697  if (max_lat < 90.0) then
698  add_np = .true.
699  jdp = jd+1
700  allocate(lat_inp(jdp))
701  lat_inp(1:jd) = lat_in(:)
702  lat_inp(jd+1) = 90.0
703  deallocate(lat_in)
704  allocate(lat_in(1:jdp))
705  lat_in(:) = lat_inp(:)
706  else
707  jdp=jd
708  endif
709  call horiz_interp_init()
710  lon_in = lon_in*pi_180
711  lat_in = lat_in*pi_180
712  allocate(x_in(id,jdp), y_in(id,jdp))
713  call meshgrid(lon_in, lat_in, x_in, y_in)
714  lon_out(:,:) = g%geoLonT(:,:)*pi_180
715  lat_out(:,:) = g%geoLatT(:,:)*pi_180
716  allocate(data_in(id,jd,kd)) ; data_in(:,:,:)=0.0
717  allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0
718  allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0
719  allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0
720  allocate(last_row(id)) ; last_row(:)=0.0
721  else
722  allocate(data_in(isd:ied,jsd:jed,kd))
723  endif
724  ! construct level cell boundaries as the mid-point between adjacent centers
725  z_edges_in(1) = 0.0
726  do k=2,kd
727  z_edges_in(k) = 0.5*(z_in(k-1)+z_in(k))
728  enddo
729  z_edges_in(kd+1) = 2.0*z_in(kd) - z_in(kd-1)
730 
731 
732  max_depth = maxval(g%bathyT)
733  call mpp_max(max_depth)
734 
735  if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
736 
737  ! roundoff = 3.0*EPSILON(missing_value)
738  roundoff = 1.e-4
739 
740  if (.not.spongedataongrid) then
741  if (is_root_pe()) &
742  call time_interp_external(fms_id, time, data_in, verbose=.true.)
743  ! loop through each data level and interpolate to model grid.
744  ! after interpolating, fill in points which will be needed
745  ! to define the layers
746  do k=1,kd
747  write(laynum,'(I8)') k ; laynum = adjustl(laynum)
748  if (is_root_pe()) then
749  tr_in(1:id,1:jd) = data_in(1:id,1:jd,k)
750  if (add_np) then
751  last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0
752  do i=1,id
753  if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then
754  pole = pole+last_row(i)
755  npole = npole+1.0
756  endif
757  enddo
758  if (npole > 0) then
759  pole=pole/npole
760  else
761  pole=missing_value
762  endif
763  tr_inp(:,1:jd) = tr_in(:,:)
764  tr_inp(:,jdp) = pole
765  else
766  tr_inp(:,:) = tr_in(:,:)
767  endif
768  endif
769 
770  call mpp_sync()
771  call mpp_broadcast(tr_inp, id*jdp, root_pe())
772  call mpp_sync_self()
773 
774  mask_in=0.0
775 
776  do j=1,jdp ; do i=1,id
777  if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then
778  mask_in(i,j)=1.0
779  tr_inp(i,j) = tr_inp(i,j) * conversion
780  else
781  tr_inp(i,j) = missing_value
782  endif
783  enddo ; enddo
784 
785  ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
786  if (k == 1) then
787  call horiz_interp_new(interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), &
788  interp_method='bilinear', src_modulo=.true.)
789  endif
790 
791  if (debug) then
792  call mystats(tr_in, missing_value, 1, id, 1, jd, k, 'Tracer from file')
793  endif
794 
795  tr_out(:,:) = 0.0
796 
797  call horiz_interp(interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, &
798  new_missing_handle=.true.)
799 
800  mask_out(:,:) = 1.0
801  do j=js,je ; do i=is,ie
802  if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j) = 0.
803  enddo ; enddo
804 
805  fill(:,:) = 0.0 ; good(:,:) = 0.0
806 
807  npoints = 0 ; varavg = 0.
808  do j=js,je ; do i=is,ie
809  if (mask_out(i,j) < 1.0) then
810  tr_out(i,j) = missing_value
811  else
812  good(i,j) = 1.0
813  npoints = npoints + 1
814  varavg = varavg + tr_out(i,j)
815  endif
816  if ((g%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= g%bathyT(i,j)) .and. &
817  (mask_out(i,j) < 1.0)) &
818  fill(i,j)=1.0
819  enddo ; enddo
820  call pass_var(fill, g%Domain)
821  call pass_var(good, g%Domain)
822 
823  if (debug) then
824  call mystats(tr_out, missing_value, is, ie, js, je, k, 'variable from horiz_interp()')
825  endif
826 
827  ! Horizontally homogenize data to produce perfectly "flat" initial conditions
828  if (PRESENT(homogenize)) then ; if (homogenize) then
829  call sum_across_pes(npoints)
830  call sum_across_pes(varavg)
831  if (npoints>0) then
832  varavg = varavg/real(npoints)
833  endif
834  tr_out(:,:) = varavg
835  endif ; endif
836 
837  ! tr_out contains input z-space data on the model grid with missing values
838  ! now fill in missing values using "ICE-nine" algorithm.
839  tr_outf(:,:) = tr_out(:,:)
840  if (k==1) tr_prev(:,:) = tr_outf(:,:)
841  good2(:,:) = good(:,:)
842  fill2(:,:) = fill(:,:)
843 
844  call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, smooth=.true.)
845 
846 ! if (debug) then
847 ! call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI)
848 ! endif
849 
850 ! call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()')
851 
852  tr_z(:,:,k) = tr_outf(:,:)*g%mask2dT(:,:)
853  mask_z(:,:,k) = good2(:,:) + fill2(:,:)
854  tr_prev(:,:) = tr_z(:,:,k)
855 
856  if (debug) then
857  call hchksum(tr_prev,'field after fill ',g%HI)
858  endif
859 
860  enddo ! kd
861  else
862  call time_interp_external(fms_id, time, data_in, verbose=.true.)
863  do k=1,kd
864  do j=js,je
865  do i=is,ie
866  tr_z(i,j,k)=data_in(i,j,k)
867  if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
868  enddo
869  enddo
870  enddo
871  endif
872 

◆ 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]filenamePath to file containing tracer to be interpolated.
[in]varnamName of tracer in filee.
[in]conversionConversion factor for tracer.
[in]recnumRecord number of tracer to be read.
[in,out]gGrid object
tr_zpointer to allocatable tracer array on local model grid and input-file vertical levels.
mask_zpointer to allocatable tracer mask array on local model grid and input-file vertical levels.
z_inCell grid values for input data.
z_edges_inCell grid edge values for input data.
[out]missing_valueThe missing value in the returned array.
[in]reentrant_xIf true, this grid is reentrant in the x-direction
[in]tripolar_nIf true, this is a northern tripolar grid
[in]homogenizeIf present and true, horizontally homogenize data to produce perfectly "flat" initial conditions
[in]m_to_zA 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.

268 
269  character(len=*), intent(in) :: filename !< Path to file containing tracer to be
270  !! interpolated.
271  character(len=*), intent(in) :: varnam !< Name of tracer in filee.
272  real, intent(in) :: conversion !< Conversion factor for tracer.
273  integer, intent(in) :: recnum !< Record number of tracer to be read.
274  type(ocean_grid_type), intent(inout) :: G !< Grid object
275  real, allocatable, dimension(:,:,:) :: tr_z !< pointer to allocatable tracer array on local
276  !! model grid and input-file vertical levels.
277  real, allocatable, dimension(:,:,:) :: mask_z !< pointer to allocatable tracer mask array on
278  !! local model grid and input-file vertical levels.
279  real, allocatable, dimension(:) :: z_in !< Cell grid values for input data.
280  real, allocatable, dimension(:) :: z_edges_in !< Cell grid edge values for input data.
281  real, intent(out) :: missing_value !< The missing value in the returned array.
282  logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction
283  logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid
284  logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data
285  !! to produce perfectly "flat" initial conditions
286  real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units
287  !! of depth. If missing, G%bathyT must be in m.
288 
289  ! Local variables
290  real, dimension(:,:), allocatable :: tr_in, tr_inp ! A 2-d array for holding input data on
291  ! native horizontal grid and extended grid
292  ! with poles.
293  real, dimension(:,:), allocatable :: mask_in ! A 2-d mask for extended input grid.
294 
295  real :: PI_180
296  integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp
297  integer :: i,j,k
298  integer, dimension(4) :: start, count, dims, dim_id
299  real, dimension(:,:), allocatable :: x_in, y_in
300  real, dimension(:), allocatable :: lon_in, lat_in
301  real, dimension(:), allocatable :: lat_inp, last_row
302  real :: max_lat, min_lat, pole, max_depth, npole
303  real :: roundoff ! The magnitude of roundoff, usually ~2e-16.
304  real :: add_offset, scale_factor
305  logical :: add_np
306  character(len=8) :: laynum
307  type(horiz_interp_type) :: Interp
308  integer :: is, ie, js, je ! compute domain indices
309  integer :: isc,iec,jsc,jec ! global compute domain indices
310  integer :: isg, ieg, jsg, jeg ! global extent
311  integer :: isd, ied, jsd, jed ! data domain indices
312  integer :: id_clock_read
313  character(len=12) :: dim_name(4)
314  logical :: debug=.false.
315  real :: npoints,varAvg
316  real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out
317  real, dimension(SZI_(G),SZJ_(G)) :: good, fill
318  real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev
319  real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2
320  real, dimension(SZI_(G),SZJ_(G)) :: nlevs
321 
322  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
323  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
324  isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
325 
326  id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=clock_loop)
327 
328 
329  if (allocated(tr_z)) deallocate(tr_z)
330  if (allocated(mask_z)) deallocate(mask_z)
331  if (allocated(z_edges_in)) deallocate(z_edges_in)
332 
333  pi_180=atan(1.0)/45.
334 
335  ! Open NetCDF file and if present, extract data and spatial coordinate information
336  ! The convention adopted here requires that the data be written in (i,j,k) ordering.
337 
338  call cpu_clock_begin(id_clock_read)
339 
340  rcode = nf90_open(filename, nf90_nowrite, ncid)
341  if (rcode /= 0) call mom_error(fatal,"error opening file "//trim(filename)//&
342  " in hinterp_extrap")
343  rcode = nf90_inq_varid(ncid, varnam, varid)
344  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(varnam)//&
345  " in file "//trim(filename)//" in hinterp_extrap")
346 
347  rcode = nf90_inquire_variable(ncid, varid, ndims=ndims, dimids=dims)
348  if (rcode /= 0) call mom_error(fatal,'error inquiring dimensions hinterp_extrap')
349  if (ndims < 3) call mom_error(fatal,"Variable "//trim(varnam)//" in file "// &
350  trim(filename)//" has too few dimensions.")
351 
352  rcode = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)
353  if (rcode /= 0) call mom_error(fatal,"error reading dimension 1 data for "// &
354  trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
355  rcode = nf90_inq_varid(ncid, dim_name(1), dim_id(1))
356  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(1))//&
357  " in file "//trim(filename)//" in hinterp_extrap")
358  rcode = nf90_inquire_dimension(ncid, dims(2), dim_name(2), len=jd)
359  if (rcode /= 0) call mom_error(fatal,"error reading dimension 2 data for "// &
360  trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
361  rcode = nf90_inq_varid(ncid, dim_name(2), dim_id(2))
362  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(2))//&
363  " in file "//trim(filename)//" in hinterp_extrap")
364  rcode = nf90_inquire_dimension(ncid, dims(3), dim_name(3), len=kd)
365  if (rcode /= 0) call mom_error(fatal,"error reading dimension 3 data for "// &
366  trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
367  rcode = nf90_inq_varid(ncid, dim_name(3), dim_id(3))
368  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(3))//&
369  " in file "//trim(filename)//" in hinterp_extrap")
370 
371  missing_value=0.0
372  rcode = nf90_get_att(ncid, varid, "_FillValue", missing_value)
373  if (rcode /= 0) call mom_error(fatal,"error finding missing value for "//&
374  trim(varnam)//" in file "// trim(filename)//" in hinterp_extrap")
375 
376  rcode = nf90_get_att(ncid, varid, "add_offset", add_offset)
377  if (rcode /= 0) add_offset = 0.0
378 
379  rcode = nf90_get_att(ncid, varid, "scale_factor", scale_factor)
380  if (rcode /= 0) scale_factor = 1.0
381 
382  if (allocated(lon_in)) deallocate(lon_in)
383  if (allocated(lat_in)) deallocate(lat_in)
384  if (allocated(z_in)) deallocate(z_in)
385  if (allocated(z_edges_in)) deallocate(z_edges_in)
386  if (allocated(tr_z)) deallocate(tr_z)
387  if (allocated(mask_z)) deallocate(mask_z)
388 
389  allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1))
390  allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))
391 
392  start = 1; count = 1; count(1) = id
393  rcode = nf90_get_var(ncid, dim_id(1), lon_in, start, count)
394  if (rcode /= 0) call mom_error(fatal,"error reading dimension 1 values for var_name "// &
395  trim(varnam)//",dim_name "//trim(dim_name(1))//" in file "// trim(filename)//" in hinterp_extrap")
396  start = 1; count = 1; count(1) = jd
397  rcode = nf90_get_var(ncid, dim_id(2), lat_in, start, count)
398  if (rcode /= 0) call mom_error(fatal,"error reading dimension 2 values for var_name "// &
399  trim(varnam)//",dim_name "//trim(dim_name(2))//" in file "// trim(filename)//" in hinterp_extrap")
400  start = 1; count = 1; count(1) = kd
401  rcode = nf90_get_var(ncid, dim_id(3), z_in, start, count)
402  if (rcode /= 0) call mom_error(fatal,"error reading dimension 3 values for var_name "// &
403  trim(varnam//",dim_name "//trim(dim_name(3)))//" in file "// trim(filename)//" in hinterp_extrap")
404 
405  call cpu_clock_end(id_clock_read)
406 
407  if (present(m_to_z)) then ; do k=1,kd ; z_in(k) = m_to_z * z_in(k) ; enddo ; endif
408 
409  ! extrapolate the input data to the north pole using the northerm-most latitude
410 
411  max_lat = maxval(lat_in)
412  add_np=.false.
413  if (max_lat < 90.0) then
414  add_np=.true.
415  jdp=jd+1
416  allocate(lat_inp(jdp))
417  lat_inp(1:jd)=lat_in(:)
418  lat_inp(jd+1)=90.0
419  deallocate(lat_in)
420  allocate(lat_in(1:jdp))
421  lat_in(:)=lat_inp(:)
422  else
423  jdp=jd
424  endif
425 
426  ! construct level cell boundaries as the mid-point between adjacent centers
427 
428  z_edges_in(1) = 0.0
429  do k=2,kd
430  z_edges_in(k)=0.5*(z_in(k-1)+z_in(k))
431  enddo
432  z_edges_in(kd+1)=2.0*z_in(kd) - z_in(kd-1)
433 
434  call horiz_interp_init()
435 
436  lon_in = lon_in*pi_180
437  lat_in = lat_in*pi_180
438  allocate(x_in(id,jdp),y_in(id,jdp))
439  call meshgrid(lon_in,lat_in, x_in, y_in)
440 
441  lon_out(:,:) = g%geoLonT(:,:)*pi_180
442  lat_out(:,:) = g%geoLatT(:,:)*pi_180
443 
444 
445  allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0
446  allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0
447  allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0
448  allocate(last_row(id)) ; last_row(:)=0.0
449 
450  max_depth = maxval(g%bathyT)
451  call mpp_max(max_depth)
452 
453  if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
454 
455  roundoff = 3.0*epsilon(missing_value)
456 
457  ! loop through each data level and interpolate to model grid.
458  ! after interpolating, fill in points which will be needed
459  ! to define the layers
460  do k=1,kd
461  write(laynum,'(I8)') k ; laynum = adjustl(laynum)
462 
463  if (is_root_pe()) then
464  start = 1 ; start(3) = k ; count(:) = 1 ; count(1) = id ; count(2) = jd
465  rcode = nf90_get_var(ncid,varid, tr_in, start, count)
466  if (rcode /= 0) call mom_error(fatal,"hinterp_and_extract_from_Fie: "//&
467  "error reading level "//trim(laynum)//" of variable "//&
468  trim(varnam)//" in file "// trim(filename))
469 
470  if (add_np) then
471  last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0
472  do i=1,id
473  if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then
474  pole = pole+last_row(i)
475  npole = npole+1.0
476  endif
477  enddo
478  if (npole > 0) then
479  pole=pole/npole
480  else
481  pole=missing_value
482  endif
483  tr_inp(:,1:jd) = tr_in(:,:)
484  tr_inp(:,jdp) = pole
485  else
486  tr_inp(:,:) = tr_in(:,:)
487  endif
488  endif
489 
490  call mpp_sync()
491  call mpp_broadcast(tr_inp, id*jdp, root_pe())
492  call mpp_sync_self()
493 
494  mask_in=0.0
495 
496  do j=1,jdp
497  do i=1,id
498  if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then
499  mask_in(i,j) = 1.0
500  tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion
501  else
502  tr_inp(i,j) = missing_value
503  endif
504  enddo
505  enddo
506 
507 ! call fms routine horiz_interp to interpolate input level data to model horizontal grid
508 
509  if (k == 1) then
510  call horiz_interp_new(interp,x_in,y_in,lon_out(is:ie,js:je),lat_out(is:ie,js:je), &
511  interp_method='bilinear',src_modulo=.true.)
512  endif
513 
514  if (debug) then
515  call mystats(tr_inp,missing_value, is,ie,js,je,k,'Tracer from file')
516  endif
517 
518  tr_out(:,:) = 0.0
519 
520  call horiz_interp(interp,tr_inp,tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.)
521 
522  mask_out=1.0
523  do j=js,je
524  do i=is,ie
525  if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j)=0.
526  enddo
527  enddo
528 
529  fill = 0.0; good = 0.0
530 
531  npoints = 0 ; varavg = 0.
532  do j=js,je
533  do i=is,ie
534  if (mask_out(i,j) < 1.0) then
535  tr_out(i,j)=missing_value
536  else
537  good(i,j)=1.0
538  npoints = npoints + 1
539  varavg = varavg + tr_out(i,j)
540  endif
541  if (g%mask2dT(i,j) == 1.0 .and. z_edges_in(k) <= g%bathyT(i,j) .and. mask_out(i,j) < 1.0) &
542  fill(i,j)=1.0
543  enddo
544  enddo
545  call pass_var(fill,g%Domain)
546  call pass_var(good,g%Domain)
547 
548  if (debug) then
549  call mystats(tr_out,missing_value, is,ie,js,je,k,'variable from horiz_interp()')
550  endif
551 
552  ! Horizontally homogenize data to produce perfectly "flat" initial conditions
553  if (PRESENT(homogenize)) then
554  if (homogenize) then
555  call sum_across_pes(npoints)
556  call sum_across_pes(varavg)
557  if (npoints>0) then
558  varavg = varavg/real(npoints)
559  endif
560  tr_out(:,:) = varavg
561  endif
562  endif
563 
564  ! tr_out contains input z-space data on the model grid with missing values
565  ! now fill in missing values using "ICE-nine" algorithm.
566  tr_outf(:,:) = tr_out(:,:)
567  if (k==1) tr_prev(:,:) = tr_outf(:,:)
568  good2(:,:) = good(:,:)
569  fill2(:,:) = fill(:,:)
570 
571  call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, smooth=.true.)
572  call mystats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()')
573 
574  tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
575  mask_z(:,:,k) = good2(:,:) + fill2(:,:)
576 
577  tr_prev(:,:) = tr_z(:,:,k)
578 
579  if (debug) then
580  call hchksum(tr_prev,'field after fill ',g%HI)
581  endif
582 
583  enddo ! kd
584 
The documentation for this interface was generated from the following file: