82 use diag_axis_mod,
only : get_diag_axis_name
83 use diag_manager_mod,
only : diag_axis_init
86 implicit none ;
private
104 logical :: configured = .false.
105 logical :: initialized = .false.
106 logical :: used = .false.
107 integer :: vertical_coord = 0
108 character(len=10) :: vertical_coord_name =
''
109 character(len=16) :: diag_coord_name =
''
110 character(len=8) :: diag_module_suffix =
''
114 real,
dimension(:,:,:),
allocatable :: h
115 real,
dimension(:),
allocatable :: dz
116 integer :: interface_axes_id = 0
117 integer :: layer_axes_id = 0
125 character(len=*),
intent(in) :: coord_tuple
128 remap_cs%diag_module_suffix = trim(
extractword(coord_tuple, 1))
129 remap_cs%diag_coord_name = trim(
extractword(coord_tuple, 2))
130 remap_cs%vertical_coord_name = trim(
extractword(coord_tuple, 3))
131 remap_cs%vertical_coord =
coordinatemode(remap_cs%vertical_coord_name)
132 remap_cs%configured = .false.
133 remap_cs%initialized = .false.
134 remap_cs%used = .false.
144 if (
allocated(remap_cs%h))
deallocate(remap_cs%h)
145 if (
allocated(remap_cs%dz))
deallocate(remap_cs%dz)
146 remap_cs%configured = .false.
147 remap_cs%initialized = .false.
148 remap_cs%used = .false.
161 if (.not. remap_cs%used)
then
173 remap_cs%used = .true.
185 integer :: nzi(4), nzl(4), k
186 character(len=200) :: inputdir, string, filename, int_varname, layer_varname
187 character(len=40) :: mod =
"MOM_diag_remap"
188 character(len=8) :: units, expected_units
189 character(len=34) :: longname, string2
191 character(len=256) :: err_msg
194 real,
allocatable,
dimension(:) :: interfaces, layers
196 call initialize_regridding(remap_cs%regrid_cs, gv, us, gv%max_depth, param_file, mod, &
197 trim(remap_cs%vertical_coord_name),
"DIAG_COORD", trim(remap_cs%diag_coord_name))
198 call set_regrid_params(remap_cs%regrid_cs, min_thickness=0., integrate_downward_for_e=.false.)
200 remap_cs%nz = get_regrid_size(remap_cs%regrid_cs)
204 longname =
'Fraction'
207 longname =
'Target Potential Density'
214 allocate(interfaces(remap_cs%nz+1))
215 allocate(layers(remap_cs%nz))
217 interfaces(:) = getcoordinateinterfaces(remap_cs%regrid_cs)
218 layers(:) = 0.5 * ( interfaces(1:remap_cs%nz) + interfaces(2:remap_cs%nz+1) )
220 remap_cs%interface_axes_id = diag_axis_init(
lowercase(trim(remap_cs%diag_coord_name))//
'_i', &
221 interfaces, trim(units),
'z', &
222 trim(longname)//
' at interface', direction=-1)
223 remap_cs%layer_axes_id = diag_axis_init(
lowercase(trim(remap_cs%diag_coord_name))//
'_l', &
224 layers, trim(units),
'z', &
225 trim(longname)//
' at cell center', direction=-1, &
226 edges=remap_cs%interface_axes_id)
229 remap_cs%configured = .true.
231 deallocate(interfaces)
240 integer,
intent(out) :: nz
241 integer,
intent(out) :: id_layer
242 integer,
intent(out) :: id_interface
245 id_layer = remap_cs%layer_axes_id
246 id_interface = remap_cs%interface_axes_id
272 real,
dimension(:, :, :),
intent(in) :: h
273 real,
dimension(:, :, :),
intent(in) :: t
274 real,
dimension(:, :, :),
intent(in) :: s
275 type(
eos_type),
pointer :: eqn_of_state
278 real,
dimension(remap_cs%nz + 1) :: zinterfaces
279 real :: h_neglect, h_neglect_edge
280 integer :: i, j, k, nz
284 if (.not. remap_cs%configured)
then
289 if (gv%Boussinesq)
then
290 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
292 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
296 if (.not. remap_cs%initialized)
then
298 call initialize_remapping(remap_cs%remap_cs,
'PPM_IH4', boundary_extrapolation=.false.)
299 allocate(remap_cs%h(g%isd:g%ied,g%jsd:g%jed, nz))
300 remap_cs%initialized = .true.
306 do j=g%jsc-1, g%jec+1 ;
do i=g%isc-1, g%iec+1
307 if (g%mask2dT(i,j)==0.)
then
308 remap_cs%h(i,j,:) = 0.
314 gv%Z_to_H*g%bathyT(i,j), sum(h(i,j,:)), &
315 zinterfaces, zscale=gv%Z_to_H)
318 gv%Z_to_H*g%bathyT(i,j), sum(h(i,j,:)), zinterfaces)
321 us%Z_to_m*g%bathyT(i,j), h(i,j,:), t(i,j,:), s(i,j,:), &
322 eqn_of_state, zinterfaces, h_neglect, h_neglect_edge)
326 call mom_error(fatal,
"diag_remap_update: SLIGHT coordinate not coded for diagnostics yet!")
330 call mom_error(fatal,
"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!")
332 remap_cs%h(i,j,:) = zinterfaces(1:nz) - zinterfaces(2:nz+1)
339 mask, missing_value, field, remapped_field)
343 real,
dimension(:,:,:),
intent(in) :: h
344 logical,
intent(in) :: staggered_in_x
345 logical,
intent(in) :: staggered_in_y
346 real,
dimension(:,:,:),
pointer :: mask
347 real,
intent(in) :: missing_value
348 real,
dimension(:,:,:),
intent(in) :: field(:,:,:)
349 real,
dimension(:,:,:),
intent(inout) :: remapped_field
351 real,
dimension(remap_cs%nz) :: h_dest
352 real,
dimension(size(h,3)) :: h_src
353 real :: h_neglect, h_neglect_edge
354 integer :: nz_src, nz_dest
357 integer :: i_lo, i_hi, j_lo, j_hi
360 call assert(remap_cs%initialized,
'diag_remap_do_remap: remap_cs not initialized.')
361 call assert(
size(field, 3) ==
size(h, 3), &
362 'diag_remap_do_remap: Remap field and thickness z-axes do not match.')
365 if (gv%Boussinesq)
then
366 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
368 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
371 nz_src =
size(field,3)
372 nz_dest = remap_cs%nz
373 remapped_field(:,:,:) = 0.
376 shift = 0;
if (g%symmetric) shift = 1
378 if (staggered_in_x .and. .not. staggered_in_y)
then
383 i_lo = i1 - shift; i_hi = i_lo + 1
384 if (
associated(mask))
then
385 if (mask(i,j,1) == 0.) cycle
387 h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
388 h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:))
390 nz_src, h_src(:), field(i1,j,:), &
391 nz_dest, h_dest(:), remapped_field(i1,j,:), &
392 h_neglect, h_neglect_edge)
395 elseif (staggered_in_y .and. .not. staggered_in_x)
then
399 j_lo = j1 - shift; j_hi = j_lo + 1
401 if (
associated(mask))
then
402 if (mask(i,j,1) == 0.) cycle
404 h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
405 h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:))
407 nz_src, h_src(:), field(i,j1,:), &
408 nz_dest, h_dest(:), remapped_field(i,j1,:), &
409 h_neglect, h_neglect_edge)
412 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y))
then
416 if (
associated(mask))
then
417 if (mask(i,j,1) == 0.) cycle
420 h_dest(:) = remap_cs%h(i,j,:)
422 nz_src, h_src(:), field(i,j,:), &
423 nz_dest, h_dest(:), remapped_field(i,j,:), &
424 h_neglect, h_neglect_edge)
428 call assert(.false.,
'diag_remap_do_remap: Unsupported axis combination')
437 real,
dimension(:,:,:),
intent(out) :: mask
439 real,
dimension(remap_cs%nz) :: h_dest
441 logical :: mask_vanished_layers
444 call assert(remap_cs%initialized,
'diag_remap_calc_hmask: remap_cs not initialized.')
447 mask_vanished_layers = (remap_cs%vertical_coord ==
coordinatemode(
'ZSTAR'))
450 do j=g%jsc-1, g%jec+1 ;
do i=g%isc-1, g%iec+1
451 if (g%mask2dT(i,j)>0.)
then
452 if (mask_vanished_layers)
then
453 h_dest(:) = remap_cs%h(i,j,:)
457 h_tot = h_tot + h_dest(k)
460 h_err = h_err + epsilon(h_tot) * h_tot
462 if (h_dest(k)<=8.*h_err)
then
478 mask, missing_value, field, reintegrated_field)
481 real,
dimension(:,:,:),
intent(in) :: h
482 logical,
intent(in) :: staggered_in_x
483 logical,
intent(in) :: staggered_in_y
484 real,
dimension(:,:,:),
pointer :: mask
485 real,
intent(in) :: missing_value
486 real,
dimension(:,:,:),
intent(in) :: field
487 real,
dimension(:,:,:),
intent(inout) :: reintegrated_field
489 real,
dimension(remap_cs%nz) :: h_dest
490 real,
dimension(size(h,3)) :: h_src
491 integer :: nz_src, nz_dest
494 integer :: i_lo, i_hi, j_lo, j_hi
497 call assert(remap_cs%initialized,
'vertically_reintegrate_diag_field: remap_cs not initialized.')
498 call assert(
size(field, 3) ==
size(h, 3), &
499 'vertically_reintegrate_diag_field: Remap field and thickness z-axes do not match.')
501 nz_src =
size(field,3)
502 nz_dest = remap_cs%nz
503 reintegrated_field(:,:,:) = 0.
506 shift = 0;
if (g%symmetric) shift = 1
508 if (staggered_in_x .and. .not. staggered_in_y)
then
513 i_lo = i1 - shift; i_hi = i_lo + 1
514 if (
associated(mask))
then
515 if (mask(i,j,1) == 0.) cycle
517 h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
518 h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:))
520 nz_dest, h_dest, 0., reintegrated_field(i1,j,:))
523 elseif (staggered_in_y .and. .not. staggered_in_x)
then
527 j_lo = j1 - shift; j_hi = j_lo + 1
529 if (
associated(mask))
then
530 if (mask(i,j,1) == 0.) cycle
532 h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
533 h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:))
535 nz_dest, h_dest, 0., reintegrated_field(i,j1,:))
538 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y))
then
542 if (
associated(mask))
then
543 if (mask(i,j,1) == 0.) cycle
546 h_dest(:) = remap_cs%h(i,j,:)
548 nz_dest, h_dest, 0., reintegrated_field(i,j,:))
552 call assert(.false.,
'vertically_reintegrate_diag_field: Q point remapping is not coded yet.')
559 mask, missing_value, field, interpolated_field)
562 real,
dimension(:,:,:),
intent(in) :: h
563 logical,
intent(in) :: staggered_in_x
564 logical,
intent(in) :: staggered_in_y
565 real,
dimension(:,:,:),
pointer :: mask
566 real,
intent(in) :: missing_value
567 real,
dimension(:,:,:),
intent(in) :: field
568 real,
dimension(:,:,:),
intent(inout) :: interpolated_field
570 real,
dimension(remap_cs%nz) :: h_dest
571 real,
dimension(size(h,3)) :: h_src
572 integer :: nz_src, nz_dest
575 integer :: i_lo, i_hi, j_lo, j_hi
578 call assert(remap_cs%initialized,
'vertically_interpolate_diag_field: remap_cs not initialized.')
579 call assert(
size(field, 3) ==
size(h, 3)+1, &
580 'vertically_interpolate_diag_field: Remap field and thickness z-axes do not match.')
582 interpolated_field(:,:,:) = 0.
585 nz_dest = remap_cs%nz
588 shift = 0;
if (g%symmetric) shift = 1
590 if (staggered_in_x .and. .not. staggered_in_y)
then
595 i_lo = i1 - shift; i_hi = i_lo + 1
596 if (
associated(mask))
then
597 if (mask(i,j,1) == 0.) cycle
599 h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
600 h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:))
602 nz_dest, h_dest, 0., interpolated_field(i1,j,:))
605 elseif (staggered_in_y .and. .not. staggered_in_x)
then
609 j_lo = j1 - shift; j_hi = j_lo + 1
611 if (
associated(mask))
then
612 if (mask(i,j,1) == 0.) cycle
614 h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
615 h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:))
617 nz_dest, h_dest, 0., interpolated_field(i,j1,:))
620 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y))
then
624 if (
associated(mask))
then
625 if (mask(i,j,1) == 0.) cycle
628 h_dest(:) = remap_cs%h(i,j,:)
630 nz_dest, h_dest, 0., interpolated_field(i,j,:))
634 call assert(.false.,
'vertically_interpolate_diag_field: Q point remapping is not coded yet.')
641 is_layer, is_extensive, &
642 missing_value, field, averaged_field, &
646 real,
dimension(:,:,:),
intent(in) :: h
647 logical,
intent(in) :: staggered_in_x
648 logical,
intent(in) :: staggered_in_y
649 logical,
intent(in) :: is_layer
650 logical,
intent(in) :: is_extensive
651 real,
intent(in) :: missing_value
652 real,
dimension(:,:,:),
intent(in) :: field
653 real,
dimension(:),
intent(inout) :: averaged_field
654 logical,
dimension(:),
intent(inout) :: averaged_mask
657 real,
dimension(G%isc:G%iec, G%jsc:G%jec, size(field,3)) :: volume, stuff
658 real,
dimension(size(field, 3)) :: vol_sum, stuff_sum
660 integer :: i, j, k, nz
669 if (staggered_in_x .and. .not. staggered_in_y)
then
675 if (is_extensive)
then
676 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
678 volume(i,j,k) = (g%US%L_to_m**2 * g%areaCu(i,j)) * g%mask2dCu(i,j)
679 stuff(i,j,k) = volume(i,j,k) * field(i1,j,k)
682 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
684 height = 0.5 * (h(i,j,k) + h(i+1,j,k))
685 volume(i,j,k) = (g%US%L_to_m**2 * g%areaCu(i,j)) &
686 * (gv%H_to_m * height) * g%mask2dCu(i,j)
687 stuff(i,j,k) = volume(i,j,k) * field(i1,j,k)
693 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
695 volume(i,j,k) = (g%US%L_to_m**2 * g%areaCu(i,j)) * g%mask2dCu(i,j)
696 stuff(i,j,k) = volume(i,j,k) * field(i1,j,k)
700 elseif (staggered_in_y .and. .not. staggered_in_x)
then
704 if (is_extensive)
then
705 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
707 volume(i,j,k) = (g%US%L_to_m**2 * g%areaCv(i,j)) * g%mask2dCv(i,j)
708 stuff(i,j,k) = volume(i,j,k) * field(i,j1,k)
711 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
713 height = 0.5 * (h(i,j,k) + h(i,j+1,k))
714 volume(i,j,k) = (g%US%L_to_m**2 * g%areaCv(i,j)) &
715 * (gv%H_to_m * height) * g%mask2dCv(i,j)
716 stuff(i,j,k) = volume(i,j,k) * field(i,j1,k)
722 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
724 volume(i,j,k) = (g%US%L_to_m**2 * g%areaCv(i,j)) * g%mask2dCv(i,j)
725 stuff(i,j,k) = volume(i,j,k) * field(i,j1,k)
729 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y))
then
733 if (is_extensive)
then
734 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
735 if (h(i,j,k) > 0.)
then
736 volume(i,j,k) = (g%US%L_to_m**2 * g%areaT(i,j)) * g%mask2dT(i,j)
737 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
744 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
745 volume(i,j,k) = (g%US%L_to_m**2 * g%areaT(i,j)) &
746 * (gv%H_to_m * h(i,j,k)) * g%mask2dT(i,j)
747 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
753 do j=g%jsc, g%jec ;
do i=g%isc, g%iec
754 volume(i,j,k) = (g%US%L_to_m**2 * g%areaT(i,j)) * g%mask2dT(i,j)
755 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
760 call assert(.false.,
'horizontally_average_diag_field: Q point averaging is not coded yet.')
768 averaged_mask(:) = .true.
770 if (vol_sum(k) > 0.)
then
771 averaged_field(k) = stuff_sum(k) / vol_sum(k)
773 averaged_field(k) = 0.
774 averaged_mask(k) = .false.