Skip to content

Commit

Permalink
Compiles with separate grid building and remapping.
Browse files Browse the repository at this point in the history
  • Loading branch information
@nichannah
nichannah committed Jul 15, 2015
1 parent bbc5330 commit 0549766
Showing 1 changed file with 90 additions and 75 deletions.
165 changes: 90 additions & 75 deletions src/framework/MOM_diag_mediator.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,21 +143,22 @@ module MOM_diag_mediator
type(remapping_CS) :: remap_cs
type(regridding_CS) :: regrid_cs
! Nominal interface locations for Z remapping
real, dimension(:), allocatable :: zi_nom
real, dimension(:), allocatable :: zi_remap
! Nominal layer locations for Z remapping
real, dimension(:), allocatable :: zl_nom
real, dimension(:), allocatable :: zl_remap
! Number of z levels used for remapping
integer :: nz_remap

! Define z star on u, v, T grids in terms of H. These are used for diagnostic remapping.
real, dimension(:,:,:), allocatable :: h_uz
real, dimension(:,:,:), allocatable :: h_vz
real, dimension(:,:,:), allocatable :: h_Tz
! Define z star on u, v, T grids, these are the interface positions
real, dimension(:,:,:), allocatable :: zi_u
real, dimension(:,:,:), allocatable :: zi_v
real, dimension(:,:,:), allocatable :: zi_T

! Keep track of which remapping is needed for diagnostic output
logical :: do_z_remapping_on_u, do_z_remapping_on_v, do_z_remapping_on_T
logical :: remapping_initialized

! Pointer to H and G for Z remapping
! Pointer to H and G for remapping
real, dimension(:,:,:), pointer :: h => null()
type(ocean_grid_type), pointer :: G => null()

Expand Down Expand Up @@ -341,74 +342,79 @@ end subroutine set_axes_info
subroutine update_diag_target_grids(G, diag_cs)
! Build target grids for diagnostic remapping

type(ocean_grid_type), intent(inout) :: G
type(diag_ctrl), intent(inout) :: diag_cs

! Arguments:
! (inout) G - ocean grid structure.
! (inout) diag_cs - structure used to regulate diagnostic output.

real, dimension(G%nk) :: h_src
real, dimension(G%iscB:G%iecB, G%jsc:G%jec, diag_cs%nz_remap) :: z_ui
real, dimension(G%jsc:G%jec, G%iscB:G%iecB, diag_cs%nz_remap) :: z_vi
real, dimension(G%jsc:G%jec, G%jsc:G%jec, diag_cs%nz_remap) :: z_Ti
real, dimension(size(diag_cs%h, 3)) :: h_src
real :: depth
integer :: nz_src, nz_dest
integer :: i, j, k

nz_dest = diag_cs%nz_remap
nz_src = size(diag_cs%h, 3)

if ((diag_cs%do_z_remapping_on_u .or. diag_cs%do_z_remapping_on_v .or. &
diag_cs%do_z_remapping_on_T) .and. (.not. diag_cs%remapping_initialized)) then
call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized')

! Initialise remapping system
call initialize_regridding(diag_cs%nz_remap, 'Z*', 'PPM_IH4', diag_cs%regrid_cs)
call initialize_remapping(G%nk, 'PPM_IH4', diag%remap_cs)
call initialize_regridding(nz_dest, 'Z*', 'PPM_IH4', diag_cs%regrid_cs)
call initialize_remapping(nz_src, 'PPM_IH4', diag_cs%remap_cs)
call setRegriddingMinimumThickness(G%Angstrom, diag_cs%regrid_cs)
call setCoordinateResolution(diag_cs%zi_remap(2:) - diag_cs%zi_remap(:diag_cs%nz_remap), diag_cs%regrid_cs)
call setCoordinateResolution(diag_cs%zi_remap(2:) - diag_cs%zi_remap(:nz_dest), diag_cs%regrid_cs)
diag_cs%remapping_initialized = .true.
endif

! Build z-star grid on u points
if (diag_cs%do_z_remapping_on_u) then
if (.not. associated(diag%z_ui)) then
allocate(diag%h_uz(G%jsc:G%jse,G%iscB,G%iecB))
if (.not. allocated(diag_cs%zi_u)) then
allocate(diag_cs%zi_u(G%iscB:G%iecB,G%jsc:G%jec,nz_dest+1))
endif
do j=G%jsc, G%jec
do i=G%iscB, G%iecB
h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:))
depth = 0.5 * (G%bathyT(i,j) + G%bathyT(i+1,j))
call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, depth, sum(h_src(:)), z_ui(i, j, :))
z_ui(i, j, :) = -z_ui(i, j, :)
diag_cs%h_uz(:) = z_ui(2:) - z_ui(:size(z_ui, 3)-1)
call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, depth, sum(h_src(:)), diag_cs%zi_u(i, j, :))
diag_cs%zi_u(i, j, :) = -diag_cs%zi_u(i, j, :)
enddo
enddo
endif

! Build z-star grid on u points
if (diag_cs%do_z_remapping_on_v) then
if (.not. associated(diag%z_ui)) then
allocate(diag%h_vz(G%jscB:G%jseB,G%isc,G%iec))
if (.not. allocated(diag_cs%zi_v)) then
allocate(diag_cs%zi_v(G%isc:G%iec,G%jscB:G%jecB,nz_dest+1))
endif

do j=G%jscB, G%jecB
do i=G%isc, G%iec
h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:))
depth = 0.5 * (G%bathyT(i, j) + G%bathyT(i, j+1))
call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, depth, total_thickness, diag%z_vi(i, j, :))
z_vi(i, j, :) = -z_vi(i, j, :)
diag_cs%h_vz(:) = z_vi(2:) - z_vi(:size(z_vi, 3)-1)
call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, depth, sum(h_src(:)), diag_cs%zi_v(i, j, :))
diag_cs%zi_v(i, j, :) = -diag_cs%zi_v(i, j, :)
enddo
enddo
endif

! Build z-star grid on T points
if (diag_cs%do_z_remapping_on_T) then
if (.not. associated(diag%z_Ti)) then
allocate(diag%h_Tz(G%jsc:G%jse,G%isc,G%iec))
if (.not. allocated(diag_cs%zi_T)) then
allocate(diag_cs%zi_T(G%isc:G%iec,G%jsc:G%jec,nz_dest+1))
endif

do j=G%jsc, G%jec
do i=G%isc, G%iec
call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, G%bathyT(i, j), sum(h(i, j, :)), diag%z_Ti(i, j, :))
z_Ti(i, j, :) = -z_Ti(i, j, :)
diag_cs%h_Tz(:) = z_Ti(2:) - z_Ti(:size(z_Ti, 3)-1)
call buildGridZstarColumn(diag_cs%regrid_cs, nz_dest, G%bathyT(i, j), sum(diag_cs%h(i, j, :)), diag_cs%zi_T(i, j, :))
diag_cs%zi_T(i, j, :) = -diag_cs%zi_T(i, j, :)
enddo
enddo
endif

end subroutine build_diag_target_grids()
end subroutine update_diag_target_grids

function check_grid_def(filename, varname)
! Do some basic checks on the vertical grid definition file, variable
Expand Down Expand Up @@ -723,8 +729,9 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field)
! (out) remapped_field - the field argument remapped to z star coordinate

real, dimension(diag_cs%nz_remap+1) :: dz
integer :: nz_src
integer :: i, j
real, dimension(diag_cs%nz_remap) :: h_dest
integer :: nz_src, nz_dest
integer :: i, j, k

call assert(size(field, 3) == size(diag_cs%h, 3), &
'Remap field and thickness z-axes do not match.')
Expand All @@ -734,72 +741,80 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field)
nz_dest = diag_cs%nz_remap

if (is_u_axes(diag%remap_axes, diag_cs)) then
call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized')

do j=diag_cs%G%jsc, diag_cs%G%jec
do i=diag_cs%G%iscB, diag_cs%G%iecB
if (associated(diag%mask3d)) then
if (diag%mask3d(i,j, 1) == 0.) cycle
endif

call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_uz, dz)
call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, &
remapped_field(:))
h_dest(:) = diag_cs%zi_u(i, j, 2:) - diag_cs%zi_u(i, j, :size(diag_cs%zi_u, 3)-1)
call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz)
call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, &
remapped_field(i, j, :))

! Lower levels of the remapped data get squashed to follow bathymetry,
! their depth does not corrospond to the nominal depth at that level
! (and the nominal depth is below the bathymetry), so remove
do k=1, nz_dest
if (diag_cs%zi_u(i, j, k) >= diag_cs%G%bathyT(i, j)) then
remapped_field(i, j, k:nz_dest) = diag_cs%missing_value
exit
endif
enddo
enddo
enddo
elseif (is_v_axes(diag%remap_axes, diag_cs)) then
call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized')

elseif (is_v_axes(diag%remap_axes, diag_cs)) then
do j=diag_cs%G%jscB, diag_cs%G%jecB
do i=diag_cs%G%isc, diag_cs%G%iec
if (associated(diag%mask3d)) then
if (diag%mask3d(i,j, 1) == 0.) cycle
endif

call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_vz, dz)
call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, &
remapped_field(:))
enddo
enddo
elseif (is_B_axes(diag%remap_axes, diag_cs)) then
call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized')

do j=diag_cs%G%jscB, diag_cs%G%jecB
do i=diag_cs%G%iscB, diag_cs%G%iecB
if (associated(diag%mask3d)) then
if (diag%mask3d(i,j, 1) == 0.) cycle
endif
call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_Bz, dz)
call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, &
remapped_field(:))
h_dest(:) = diag_cs%zi_v(i, j, 2:) - diag_cs%zi_v(i, j, :size(diag_cs%zi_v, 3)-1)
call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz)
call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, &
remapped_field(i, j, :))
do k=1, nz_dest
if (diag_cs%zi_v(i, j, k) >= diag_cs%G%bathyT(i, j)) then
remapped_field(i, j, k:nz_dest) = diag_cs%missing_value
exit
endif
enddo
enddo
enddo
!elseif (is_B_axes(diag%remap_axes, diag_cs)) then
! do j=diag_cs%G%jscB, diag_cs%G%jecB
! do i=diag_cs%G%iscB, diag_cs%G%iecB
! if (associated(diag%mask3d)) then
! if (diag%mask3d(i,j, 1) == 0.) cycle
! endif
! h_dest(:) = diag_cs%zi_B(i, j, 2:) - diag_cs%zi_B(i, j, :size(diag_cs%zi_B, 3)-1)
! call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz)
! call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, &
! remapped_field(i, j, :))
! enddo
! enddo
else
call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized')

do j=diag_cs%G%jsc, diag_cs%G%jec
do i=diag_cs%G%isc, diag_cs%G%iec
if (associated(diag%mask3d)) then
if (diag%mask3d(i,j, 1) == 0.) cycle
endif
h_dest(:) = zi_remap(2:) - zi_remap(:size(zi_remap)-1)
call dzFromH1H2(nz_src, h(:), nz_dest, diag_cs%h_Tz, dz)
call remapping_core(remap_cs, nz_src, h(:), field(:), nz_dest, dz, &
remapped_field(:))
h_dest(:) = diag_cs%zi_T(i, j, 2:) - diag_cs%zi_T(i, j, :size(diag_cs%zi_T, 3)-1)
call dzFromH1H2(nz_src, diag_cs%h(i, j, :), nz_dest, h_dest(:), dz)
call remapping_core(diag_cs%remap_cs, nz_src, diag_cs%h(i, j, :), field(i, j, :), nz_dest, dz, &
remapped_field(i, j, :))
do k=1, nz_dest
if (diag_cs%zi_T(i, j, k) >= diag_cs%G%bathyT(i, j)) then
remapped_field(i, j, k:nz_dest) = diag_cs%missing_value
exit
endif
enddo
enddo
enddo
endif

! Lower levels of the remapped data get squashed to follow bathymetry,
! their depth does not corrospond to the nominal depth at that level
! (and the nominal depth is below the bathymetry), so remove
do k=1, nz_dest
if (zi_remap(k) >= depth) then
remapped_field(k:nz_dest) = missing_value
exit
endif
enddo

end subroutine remap_diag_to_z

subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask)
Expand Down Expand Up @@ -1452,9 +1467,9 @@ subroutine diag_mediator_init(G, param_file, diag_cs, err_msg)
! Keep a pointer to the grid, this is needed for regridding
diag_cs%G => G
diag_cs%nz_remap = -1
do_z_remapping_on_u = .false.
do_z_remapping_on_v = .false.
do_z_remapping_on_T = .false.
diag_cs%do_z_remapping_on_u = .false.
diag_cs%do_z_remapping_on_v = .false.
diag_cs%do_z_remapping_on_T = .false.

diag_cs%is = G%isc - (G%isd-1) ; diag_cs%ie = G%iec - (G%isd-1)
diag_cs%js = G%jsc - (G%jsd-1) ; diag_cs%je = G%jec - (G%jsd-1)
Expand Down

0 comments on commit 0549766

Please sign in to comment.