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 53 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, answers_2018, ongrid)
 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, answers_2018)
 Extrapolate and interpolate using a FMS time interpolation handle. More...
 

Detailed Description

Extrapolate and interpolate data.

Definition at line 53 of file MOM_horizontal_regridding.F90.

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,
logical, intent(in), optional  answers_2018 
)
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.
[in]answers_2018If true, use expressions that give the same answers as the code did in late 2018. Otherwise add parentheses for rotational symmetry.

Definition at line 633 of file MOM_horizontal_regridding.F90.

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

◆ 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,
logical, intent(in), optional  answers_2018,
logical, intent(in), optional  ongrid 
)
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.
[in]answers_2018If true, use expressions that give the same answers as the code did in late 2018. Otherwise add parentheses for rotational symmetry.
[in]ongridIf true, then data are assumed to have been interpolated to the model horizontal grid. In this case, only extrapolation is performed by this routine

Definition at line 274 of file MOM_horizontal_regridding.F90.

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