7 use mom_coms,
only : max_across_pes, min_across_pes
8 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
24 use mom_time_manager,
only : get_external_field_axes, get_external_field_missing
26 use mpp_io_mod,
only : axistype
27 use mpp_domains_mod,
only : mpp_global_field, mpp_get_compute_domain
28 use mpp_mod,
only : mpp_broadcast,mpp_root_pe,mpp_sync,mpp_sync_self
29 use mpp_mod,
only : mpp_max
30 use horiz_interp_mod,
only : horiz_interp_new, horiz_interp,horiz_interp_type
31 use horiz_interp_mod,
only : horiz_interp_init, horiz_interp_del
33 use mpp_io_mod,
only : mpp_get_axis_data
34 use mpp_io_mod,
only : mpp_single
37 implicit none ;
private
39 #include <MOM_memory.h>
60 subroutine mystats(array, missing, is, ie, js, je, k, mesg)
61 real,
dimension(:,:),
intent(in) :: array
62 real,
intent(in) :: missing
65 integer :: is,ie,js,je
68 character(len=*) :: mesg
73 character(len=120) :: lmesg
74 mina = 9.e24 ; maxa = -9.e24 ; found = .false.
76 do j=js,je ;
do i=is,ie
77 if (array(i,j) /= array(i,j)) stop
'Nan!'
78 if (abs(array(i,j)-missing) > 1.e-6*abs(missing))
then
80 mina = min(mina, array(i,j))
81 maxa = max(maxa, array(i,j))
89 call min_across_pes(mina)
90 call max_across_pes(maxa)
91 if (is_root_pe())
then
92 write(lmesg(1:120),
'(2(a,es12.4),a,i3,x,a)') &
93 'init_from_Z: min=',mina,
' max=',maxa,
' Level=',k,trim(mesg)
94 call mom_mesg(lmesg,2)
103 subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, keep_bug, debug)
107 real,
dimension(SZI_(G),SZJ_(G)), &
108 intent(inout) :: aout
109 real,
dimension(SZI_(G),SZJ_(G)), &
112 real,
dimension(SZI_(G),SZJ_(G)), &
115 real,
dimension(SZI_(G),SZJ_(G)), &
116 optional,
intent(in) :: prev
117 logical,
optional,
intent(in) :: smooth
119 integer,
optional,
intent(in) :: num_pass
120 real,
optional,
intent(in) :: relc
121 real,
optional,
intent(in) :: crit
122 logical,
optional,
intent(in) :: keep_bug
124 logical,
optional,
intent(in) :: debug
126 real,
dimension(SZI_(G),SZJ_(G)) :: b,r
127 real,
dimension(SZI_(G),SZJ_(G)) :: fill_pts, good_, good_new
129 character(len=256) :: mesg
131 real :: east,west,north,south,sor
132 real :: ge,gw,gn,gs,ngood
133 logical :: do_smooth,siena_bug
134 real :: nfill, nfill_prev
135 integer,
parameter :: num_pass_default = 10000
136 real,
parameter :: relc_default = 0.25, crit_default = 1.e-3
139 integer :: is, ie, js, je
140 real :: relax_coeff, acrit, ares
144 if (
PRESENT(debug)) debug_it=debug
146 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
148 npass = num_pass_default
149 if (
PRESENT(num_pass)) npass = num_pass
151 relax_coeff = relc_default
152 if (
PRESENT(relc)) relax_coeff = relc
155 if (
PRESENT(crit)) acrit = crit
158 if (
PRESENT(keep_bug)) siena_bug = keep_bug
161 if (
PRESENT(smooth)) do_smooth=smooth
163 fill_pts(:,:) = fill(:,:)
165 nfill = sum(fill(is:ie,js:je))
166 call sum_across_pes(nfill)
169 good_(:,:) = good(:,:)
172 do while (nfill > 0.0)
178 good_new(:,:)=good_(:,:)
180 do j=js,je ;
do i=is,ie
182 if (good_(i,j) == 1.0 .or. fill(i,j) == 0.) cycle
184 ge=good_(i+1,j) ; gw=good_(i-1,j)
185 gn=good_(i,j+1) ; gs=good_(i,j-1)
186 east=0.0 ; west=0.0 ; north=0.0 ; south=0.0
187 if (ge == 1.0) east = aout(i+1,j)*ge
188 if (gw == 1.0) west = aout(i-1,j)*gw
189 if (gn == 1.0) north = aout(i,j+1)*gn
190 if (gs == 1.0) south = aout(i,j-1)*gs
194 b(i,j)=(east+west+north+south)/ngood
202 aout(is:ie,js:je) = b(is:ie,js:je)
203 good_(is:ie,js:je) = good_new(is:ie,js:je)
205 nfill = sum(fill_pts(is:ie,js:je))
206 call sum_across_pes(nfill)
208 if (nfill == nfill_prev .and.
PRESENT(prev))
then
209 do j=js,je ;
do i=is,ie ;
if (fill_pts(i,j) == 1.0)
then
210 aout(i,j) = prev(i,j)
212 endif ;
enddo ;
enddo
213 elseif (nfill == nfill_prev)
then
214 call mom_error(warning, &
215 'Unable to fill missing points using either data at the same vertical level from a connected basin'//&
216 'or using a point from a previous vertical level. Make sure that the original data has some valid'//&
217 'data in all basins.', .true.)
218 write(mesg,*)
'nfill=',nfill
219 call mom_error(warning, mesg, .true.)
222 nfill = sum(fill_pts(is:ie,js:je))
223 call sum_across_pes(nfill)
227 if (do_smooth)
then ;
do k=1,npass
229 do j=js,je ;
do i=is,ie
230 if (fill(i,j) == 1)
then
231 east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j))
232 north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1))
233 r(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + &
234 west*aout(i-1,j)+east*aout(i+1,j) - &
235 (south+north+west+east)*aout(i,j))
245 do j=js,je ;
do i=is,ie
246 aout(i,j) = r(i,j) + aout(i,j)
247 ares = max(ares, abs(r(i,j)))
249 call max_across_pes(ares)
250 if (ares <= acrit)
exit
253 do j=js,je ;
do i=is,ie
254 if (good_(i,j) == 0.0 .and. fill_pts(i,j) == 1.0)
then
255 write(mesg,*)
'In fill_miss, fill, good,i,j= ',fill_pts(i,j),good_(i,j),i,j
256 call mom_error(warning, mesg, .true.)
257 call mom_error(fatal,
"MOM_initialize: "// &
258 "fill is true and good is false after fill_miss, how did this happen? ")
266 mask_z, z_in, z_edges_in, missing_value, reentrant_x, &
267 tripolar_n, homogenize, m_to_Z)
269 character(len=*),
intent(in) :: filename
271 character(len=*),
intent(in) :: varnam
272 real,
intent(in) :: conversion
273 integer,
intent(in) :: recnum
275 real,
allocatable,
dimension(:,:,:) :: tr_z
277 real,
allocatable,
dimension(:,:,:) :: mask_z
279 real,
allocatable,
dimension(:) :: z_in
280 real,
allocatable,
dimension(:) :: z_edges_in
281 real,
intent(out) :: missing_value
282 logical,
intent(in) :: reentrant_x
283 logical,
intent(in) :: tripolar_n
284 logical,
optional,
intent(in) :: homogenize
286 real,
optional,
intent(in) :: m_to_Z
290 real,
dimension(:,:),
allocatable :: tr_in, tr_inp
293 real,
dimension(:,:),
allocatable :: mask_in
296 integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp
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
304 real :: add_offset, scale_factor
306 character(len=8) :: laynum
307 type(horiz_interp_type) :: Interp
308 integer :: is, ie, js, je
309 integer :: isc,iec,jsc,jec
310 integer :: isg, ieg, jsg, jeg
311 integer :: isd, ied, jsd, jed
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
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
326 id_clock_read = cpu_clock_id(
'(Initialize tracer from Z) read', grain=clock_loop)
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)
338 call cpu_clock_begin(id_clock_read)
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")
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.")
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")
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")
376 rcode = nf90_get_att(ncid, varid,
"add_offset", add_offset)
377 if (rcode /= 0) add_offset = 0.0
379 rcode = nf90_get_att(ncid, varid,
"scale_factor", scale_factor)
380 if (rcode /= 0) scale_factor = 1.0
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)
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))
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")
405 call cpu_clock_end(id_clock_read)
407 if (
present(m_to_z))
then ;
do k=1,kd ; z_in(k) = m_to_z * z_in(k) ;
enddo ;
endif
411 max_lat = maxval(lat_in)
413 if (max_lat < 90.0)
then
416 allocate(lat_inp(jdp))
417 lat_inp(1:jd)=lat_in(:)
420 allocate(lat_in(1:jdp))
430 z_edges_in(k)=0.5*(z_in(k-1)+z_in(k))
432 z_edges_in(kd+1)=2.0*z_in(kd) - z_in(kd-1)
434 call horiz_interp_init()
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)
441 lon_out(:,:) = g%geoLonT(:,:)*pi_180
442 lat_out(:,:) = g%geoLatT(:,:)*pi_180
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
450 max_depth = maxval(g%bathyT)
451 call mpp_max(max_depth)
453 if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
455 roundoff = 3.0*epsilon(missing_value)
461 write(laynum,
'(I8)') k ; laynum = adjustl(laynum)
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))
471 last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0
473 if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value))
then
474 pole = pole+last_row(i)
483 tr_inp(:,1:jd) = tr_in(:,:)
486 tr_inp(:,:) = tr_in(:,:)
491 call mpp_broadcast(tr_inp, id*jdp, root_pe())
498 if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value))
then
500 tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion
502 tr_inp(i,j) = missing_value
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.)
515 call mystats(tr_inp,missing_value, is,ie,js,je,k,
'Tracer from file')
520 call horiz_interp(interp,tr_inp,tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.)
525 if (abs(tr_out(i,j)-missing_value) < abs(roundoff*missing_value)) mask_out(i,j)=0.
529 fill = 0.0; good = 0.0
531 npoints = 0 ; varavg = 0.
534 if (mask_out(i,j) < 1.0)
then
535 tr_out(i,j)=missing_value
538 npoints = npoints + 1
539 varavg = varavg + tr_out(i,j)
541 if (g%mask2dT(i,j) == 1.0 .and. z_edges_in(k) <= g%bathyT(i,j) .and. mask_out(i,j) < 1.0) &
549 call mystats(tr_out,missing_value, is,ie,js,je,k,
'variable from horiz_interp()')
553 if (
PRESENT(homogenize))
then
555 call sum_across_pes(npoints)
556 call sum_across_pes(varavg)
558 varavg = varavg/real(npoints)
566 tr_outf(:,:) = tr_out(:,:)
567 if (k==1) tr_prev(:,:) = tr_outf(:,:)
568 good2(:,:) = good(:,:)
569 fill2(:,:) = fill(:,:)
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()')
574 tr_z(:,:,k) = tr_outf(:,:) * g%mask2dT(:,:)
575 mask_z(:,:,k) = good2(:,:) + fill2(:,:)
577 tr_prev(:,:) = tr_z(:,:,k)
580 call hchksum(tr_prev,
'field after fill ',g%HI)
589 z_in, z_edges_in, missing_value, reentrant_x, &
590 tripolar_n, homogenize, spongeOngrid, m_to_Z)
592 integer,
intent(in) :: fms_id
593 type(time_type),
intent(in) :: Time
594 real,
intent(in) :: conversion
596 real,
allocatable,
dimension(:,:,:) :: tr_z
598 real,
allocatable,
dimension(:,:,:) :: mask_z
600 real,
allocatable,
dimension(:) :: z_in
601 real,
allocatable,
dimension(:) :: z_edges_in
602 real,
intent(out) :: missing_value
603 logical,
intent(in) :: reentrant_x
604 logical,
intent(in) :: tripolar_n
605 logical,
optional,
intent(in) :: homogenize
607 logical,
optional,
intent(in) :: spongeOngrid
608 real,
optional,
intent(in) :: m_to_Z
612 real,
dimension(:,:),
allocatable :: tr_in,tr_inp
615 real,
dimension(:,:,:),
allocatable :: data_in
617 real,
dimension(:,:),
allocatable :: mask_in
620 integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp
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
629 character(len=8) :: laynum
630 type(horiz_interp_type) :: Interp
631 type(axistype),
dimension(4) :: axes_data
632 integer :: is, ie, js, je
633 integer :: isc,iec,jsc,jec
634 integer :: isg, ieg, jsg, jeg
635 integer :: isd, ied, jsd, jed
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
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
652 id_clock_read = cpu_clock_id(
'(Initialize tracer from Z) read', grain=clock_loop)
659 call cpu_clock_begin(id_clock_read)
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)
669 axes_data = get_external_field_axes(fms_id)
671 id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3)
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)
681 allocate(z_in(kd),z_edges_in(kd+1))
683 allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))
685 call mpp_get_axis_data(axes_data(3), z_in)
687 if (
present(m_to_z))
then ;
do k=1,kd ; z_in(k) = m_to_z * z_in(k) ;
enddo ;
endif
689 call cpu_clock_end(id_clock_read)
691 missing_value = get_external_field_missing(fms_id)
693 if (.not. spongedataongrid)
then
695 max_lat = maxval(lat_in)
697 if (max_lat < 90.0)
then
700 allocate(lat_inp(jdp))
701 lat_inp(1:jd) = lat_in(:)
704 allocate(lat_in(1:jdp))
705 lat_in(:) = lat_inp(:)
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
722 allocate(data_in(isd:ied,jsd:jed,kd))
727 z_edges_in(k) = 0.5*(z_in(k-1)+z_in(k))
729 z_edges_in(kd+1) = 2.0*z_in(kd) - z_in(kd-1)
732 max_depth = maxval(g%bathyT)
733 call mpp_max(max_depth)
735 if (z_edges_in(kd+1)<max_depth) z_edges_in(kd+1)=max_depth
740 if (.not.spongedataongrid)
then
742 call time_interp_external(fms_id, time, data_in, verbose=.true.)
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)
751 last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0
753 if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value))
then
754 pole = pole+last_row(i)
763 tr_inp(:,1:jd) = tr_in(:,:)
766 tr_inp(:,:) = tr_in(:,:)
771 call mpp_broadcast(tr_inp, id*jdp, root_pe())
776 do j=1,jdp ;
do i=1,id
777 if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value))
then
779 tr_inp(i,j) = tr_inp(i,j) * conversion
781 tr_inp(i,j) = missing_value
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.)
792 call mystats(tr_in, missing_value, 1, id, 1, jd, k,
'Tracer from file')
797 call horiz_interp(interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, &
798 new_missing_handle=.true.)
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.
805 fill(:,:) = 0.0 ; good(:,:) = 0.0
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
813 npoints = npoints + 1
814 varavg = varavg + tr_out(i,j)
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)) &
824 call mystats(tr_out, missing_value, is, ie, js, je, k,
'variable from horiz_interp()')
828 if (
PRESENT(homogenize))
then ;
if (homogenize)
then
829 call sum_across_pes(npoints)
830 call sum_across_pes(varavg)
832 varavg = varavg/real(npoints)
839 tr_outf(:,:) = tr_out(:,:)
840 if (k==1) tr_prev(:,:) = tr_outf(:,:)
841 good2(:,:) = good(:,:)
842 fill2(:,:) = fill(:,:)
844 call fill_miss_2d(tr_outf, good2, fill2, tr_prev, g, smooth=.true.)
852 tr_z(:,:,k) = tr_outf(:,:)*g%mask2dT(:,:)
853 mask_z(:,:,k) = good2(:,:) + fill2(:,:)
854 tr_prev(:,:) = tr_z(:,:,k)
857 call hchksum(tr_prev,
'field after fill ',g%HI)
862 call time_interp_external(fms_id, time, data_in, verbose=.true.)
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.
877 real,
dimension(:),
intent(in) :: x
878 real,
dimension(:),
intent(in) :: y
879 real,
dimension(size(x,1),size(y,1)),
intent(inout) :: x_T
880 real,
dimension(size(x,1),size(y,1)),
intent(inout) :: y_T
884 ni=
size(x,1) ; nj=
size(y,1)
886 do j=1,nj ;
do i=1,ni
890 do j=1,nj ;
do i=1,ni
900 integer,
dimension(:,:),
intent(in) :: m
901 logical,
intent(in) :: cyclic_x
902 logical,
intent(in) :: tripolar_n
903 integer,
dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp
905 real,
dimension(size(m,1),size(m,2)) :: m_real
906 real,
dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real
918 real,
dimension(:,:),
intent(in) :: m
919 logical,
intent(in) :: cyclic_x
920 logical,
intent(in) :: tripolar_n
921 real,
dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp
925 ni=
size(m,1); nj=
size(m,2)
930 mp(0,1:nj)=m(ni,1:nj)
931 mp(ni+1,1:nj)=m(1,1:nj)
934 mp(ni+1,1:nj)=m(ni,1:nj)
940 mp(i,nj+1)=m(ni-i+1,nj)
943 mp(1:ni,nj+1)=m(1:ni,nj)
954 subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n)
955 real,
dimension(:,:),
intent(inout) :: zi
956 integer,
dimension(size(zi,1),size(zi,2)),
intent(in) :: fill
957 integer,
dimension(size(zi,1),size(zi,2)),
intent(in) :: bad
958 real,
intent(in) :: sor
959 integer,
intent(in) :: niter
960 logical,
intent(in) :: cyclic_x
961 logical,
intent(in) :: tripolar_n
964 real,
dimension(size(zi,1),size(zi,2)) :: res, m
965 integer,
dimension(size(zi,1),size(zi,2),4) :: B
966 real,
dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp
967 integer,
dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm
972 ni=
size(zi,1) ; nj=
size(zi,2)
982 if (fill(i,j) == 1)
then
983 b(i,j,1)=1-nm(i+1,j);b(i,j,2)=1-nm(i-1,j)
984 b(i,j,3)=1-nm(i,j+1);b(i,j,4)=1-nm(i,j-1)
992 if (fill(i,j) == 1)
then
993 bsum = real(b(i,j,1)+b(i,j,2)+b(i,j,3)+b(i,j,4))
995 res(i,j)=isum*(b(i,j,1)*mp(i+1,j)+b(i,j,2)*mp(i-1,j)+&
996 b(i,j,3)*mp(i,j+1)+b(i,j,4)*mp(i,j-1)) - mp(i,j)
1000 res(:,:)=res(:,:)*sor
1004 mp(i,j)=mp(i,j)+res(i,j)
1008 zi(:,:)=mp(1:ni,1:nj)