Extrapolate and interpolate from a file record.
271 character(len=*),
intent(in) :: filename
273 character(len=*),
intent(in) :: varnam
274 real,
intent(in) :: conversion
275 integer,
intent(in) :: recnum
276 type(ocean_grid_type),
intent(inout) :: G
277 real,
allocatable,
dimension(:,:,:) :: tr_z
279 real,
allocatable,
dimension(:,:,:) :: mask_z
281 real,
allocatable,
dimension(:) :: z_in
282 real,
allocatable,
dimension(:) :: z_edges_in
283 real,
intent(out) :: missing_value
284 logical,
intent(in) :: reentrant_x
285 logical,
intent(in) :: tripolar_n
286 logical,
optional,
intent(in) :: homogenize
288 real,
optional,
intent(in) :: m_to_Z
292 real,
dimension(:,:),
allocatable :: tr_in, tr_inp
295 real,
dimension(:,:),
allocatable :: mask_in
298 integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp
300 integer,
dimension(4) :: start, count, dims, dim_id
301 real,
dimension(:,:),
allocatable :: x_in, y_in
302 real,
dimension(:),
allocatable :: lon_in, lat_in
303 real,
dimension(:),
allocatable :: lat_inp, last_row
304 real :: max_lat, min_lat, pole, max_depth, npole
306 real :: add_offset, scale_factor
308 character(len=8) :: laynum
309 type(horiz_interp_type) :: Interp
310 integer :: is, ie, js, je
311 integer :: isc,iec,jsc,jec
312 integer :: isg, ieg, jsg, jeg
313 integer :: isd, ied, jsd, jed
314 integer :: id_clock_read
315 character(len=12) :: dim_name(4)
316 logical :: debug=.false.
317 real :: npoints,varAvg
318 real,
dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out
319 real,
dimension(SZI_(G),SZJ_(G)) :: good, fill
320 real,
dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev
321 real,
dimension(SZI_(G),SZJ_(G)) :: good2,fill2
322 real,
dimension(SZI_(G),SZJ_(G)) :: nlevs
324 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
325 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
326 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
328 id_clock_read = cpu_clock_id(
'(Initialize tracer from Z) read', grain=clock_loop)
331 if (
allocated(tr_z))
deallocate(tr_z)
332 if (
allocated(mask_z))
deallocate(mask_z)
333 if (
allocated(z_edges_in))
deallocate(z_edges_in)
340 call cpu_clock_begin(id_clock_read)
343 rcode = nf90_open(filename, nf90_nowrite, ncid)
344 if (rcode /= 0)
call mom_error(fatal,
"error opening file "//trim(filename)//&
345 " in hinterp_extrap")
346 rcode = nf90_inq_varid(ncid, varnam, varid)
347 if (rcode /= 0)
call mom_error(fatal,
"error finding variable "//trim(varnam)//&
348 " in file "//trim(filename)//
" in hinterp_extrap")
350 rcode = nf90_inquire_variable(ncid, varid, ndims=ndims, dimids=dims)
351 if (rcode /= 0)
call mom_error(fatal,
'error inquiring dimensions hinterp_extrap')
352 if (ndims < 3)
call mom_error(fatal,
"Variable "//trim(varnam)//
" in file "// &
353 trim(filename)//
" has too few dimensions.")
355 rcode = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)
356 if (rcode /= 0)
call mom_error(fatal,
"error reading dimension 1 data for "// &
357 trim(varnam)//
" in file "// trim(filename)//
" in hinterp_extrap")
358 rcode = nf90_inq_varid(ncid, dim_name(1), dim_id(1))
359 if (rcode /= 0)
call mom_error(fatal,
"error finding variable "//trim(dim_name(1))//&
360 " in file "//trim(filename)//
" in hinterp_extrap")
361 rcode = nf90_inquire_dimension(ncid, dims(2), dim_name(2), len=jd)
362 if (rcode /= 0)
call mom_error(fatal,
"error reading dimension 2 data for "// &
363 trim(varnam)//
" in file "// trim(filename)//
" in hinterp_extrap")
364 rcode = nf90_inq_varid(ncid, dim_name(2), dim_id(2))
365 if (rcode /= 0)
call mom_error(fatal,
"error finding variable "//trim(dim_name(2))//&
366 " in file "//trim(filename)//
" in hinterp_extrap")
367 rcode = nf90_inquire_dimension(ncid, dims(3), dim_name(3), len=kd)
368 if (rcode /= 0)
call mom_error(fatal,
"error reading dimension 3 data for "// &
369 trim(varnam)//
" in file "// trim(filename)//
" in hinterp_extrap")
370 rcode = nf90_inq_varid(ncid, dim_name(3), dim_id(3))
371 if (rcode /= 0)
call mom_error(fatal,
"error finding variable "//trim(dim_name(3))//&
372 " in file "//trim(filename)//
" in hinterp_extrap")
376 rcode = nf90_get_att(ncid, varid,
"_FillValue", missing_value)
377 if (rcode /= 0)
call mom_error(fatal,
"error finding missing value for "//&
378 trim(varnam)//
" in file "// trim(filename)//
" in hinterp_extrap")
380 rcode = nf90_get_att(ncid, varid,
"add_offset", add_offset)
381 if (rcode /= 0) add_offset = 0.0
383 rcode = nf90_get_att(ncid, varid,
"scale_factor", scale_factor)
384 if (rcode /= 0) scale_factor = 1.0
387 if (
allocated(lon_in))
deallocate(lon_in)
388 if (
allocated(lat_in))
deallocate(lat_in)
389 if (
allocated(z_in))
deallocate(z_in)
390 if (
allocated(z_edges_in))
deallocate(z_edges_in)
391 if (
allocated(tr_z))
deallocate(tr_z)
392 if (
allocated(mask_z))
deallocate(mask_z)
395 allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1))
396 allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))
398 start = 1; count = 1; count(1) = id
399 rcode = nf90_get_var(ncid, dim_id(1), lon_in, start, count)
400 if (rcode /= 0)
call mom_error(fatal,
"error reading dimension 1 values for var_name "// &
401 trim(varnam)//
",dim_name "//trim(dim_name(1))//
" in file "// trim(filename)//
" in hinterp_extrap")
402 start = 1; count = 1; count(1) = jd
403 rcode = nf90_get_var(ncid, dim_id(2), lat_in, start, count)
404 if (rcode /= 0)
call mom_error(fatal,
"error reading dimension 2 values for var_name "// &
405 trim(varnam)//
",dim_name "//trim(dim_name(2))//
" in file "// trim(filename)//
" in hinterp_extrap")
406 start = 1; count = 1; count(1) = kd
407 rcode = nf90_get_var(ncid, dim_id(3), z_in, start, count)
408 if (rcode /= 0)
call mom_error(fatal,
"error reading dimension 3 values for var_name "// &
409 trim(varnam//
",dim_name "//trim(dim_name(3)))//
" in file "// trim(filename)//
" in hinterp_extrap")
411 call cpu_clock_end(id_clock_read)
413 if (
present(m_to_z))
then ;
do k=1,kd ; z_in(k) = m_to_z * z_in(k) ;
enddo ;
endif
417 max_lat = maxval(lat_in)
419 if (max_lat < 90.0)
then
422 allocate(lat_inp(jdp))
423 lat_inp(1:jd)=lat_in(:)
426 allocate(lat_in(1:jdp))
436 z_edges_in(k)=0.5*(z_in(k-1)+z_in(k))
438 z_edges_in(kd+1)=2.0*z_in(kd) - z_in(kd-1)
440 call horiz_interp_init()
442 lon_in = lon_in*pi_180
443 lat_in = lat_in*pi_180
444 allocate(x_in(id,jdp),y_in(id,jdp))
445 call meshgrid(lon_in,lat_in, x_in, y_in)
447 lon_out(:,:) = g%geoLonT(:,:)*pi_180
448 lat_out(:,:) = g%geoLatT(:,:)*pi_180
451 allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0
452 allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0
453 allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0
454 allocate(last_row(id)) ; last_row(:)=0.0
456 max_depth = maxval(g%bathyT)
457 call mpp_max(max_depth)
459 if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
466 roundoff = 3.0*epsilon(missing_value)
470 write(laynum,
'(I8)') k ; laynum = adjustl(laynum)
472 if (is_root_pe())
then
473 start = 1 ; start(3) = k ; count(:) = 1 ; count(1) = id ; count(2) = jd
474 rcode = nf90_get_var(ncid,varid, tr_in, start, count)
475 if (rcode /= 0)
call mom_error(fatal,
"hinterp_and_extract_from_Fie: "//&
476 "error reading level "//trim(laynum)//
" of variable "//&
477 trim(varnam)//
" in file "// trim(filename))
480 last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0
482 if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value))
then
483 pole = pole+last_row(i)
492 tr_inp(:,1:jd) = tr_in(:,:)
495 tr_inp(:,:) = tr_in(:,:)
501 call mpp_broadcast(tr_inp, id*jdp, root_pe())
508 if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value))
then
510 tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion
512 tr_inp(i,j) = missing_value
522 call horiz_interp_new(interp,x_in,y_in,lon_out(is:ie,js:je),lat_out(is:ie,js:je), &
523 interp_method=
'bilinear',src_modulo=reentrant_x)
527 call mystats(tr_inp,missing_value, is,ie,js,je,k,
'Tracer from file')
532 call horiz_interp(interp,tr_inp,tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.)
537 if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j)=0.
541 fill = 0.0; good = 0.0
543 npoints = 0 ; varavg = 0.
546 if (mask_out(i,j) < 1.0)
then
547 tr_out(i,j)=missing_value
550 npoints = npoints + 1
551 varavg = varavg + tr_out(i,j)
553 if (g%mask2dT(i,j) == 1.0 .and. z_edges_in(k) <= g%bathyT(i,j) .and. mask_out(i,j) < 1.0) &
557 call pass_var(fill,g%Domain)
558 call pass_var(good,g%Domain)
561 call mystats(tr_out,missing_value, is,ie,js,je,k,
'variable from horiz_interp()')
565 if (
PRESENT(homogenize))
then
567 call sum_across_pes(npoints)
568 call sum_across_pes(varavg)
570 varavg = varavg/real(npoints)
579 tr_outf(:,:) = tr_out(:,:)
580 if (k==1) tr_prev(:,:) = tr_outf(:,:)
581 good2(:,:) = good(:,:)
582 fill2(:,:) = fill(:,:)
584 call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, smooth=.true.)
585 call mystats(tr_outf, missing_value, is, ie, js, je, k,
'field from fill_miss_2d()')
587 tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
588 mask_z(:,:,k) = good2(:,:) + fill2(:,:)
590 tr_prev(:,:) = tr_z(:,:,k)
593 call hchksum(tr_prev,
'field after fill ',g%HI)