From 781e94d10b0bdf09d08bc2a41ab6837a2e8c8551 Mon Sep 17 00:00:00 2001 From: Nicholas Hannah Date: Thu, 7 Jul 2016 11:15:58 +1000 Subject: [PATCH] Prototyping diagnostic regridding code. #334 --- src/framework/MOM_diag_mediator.F90 | 341 ++++++++++++++++++++++------ 1 file changed, 276 insertions(+), 65 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 8cc208a232..0c1d6bc686 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -42,7 +42,10 @@ module MOM_diag_mediator use MOM_remapping, only : remapping_CS, initialize_remapping, dzFromH1H2 use MOM_remapping, only : remapping_core_h use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution -use MOM_regridding, only : build_zstar_column, set_regrid_params +use MOM_regridding, only : build_zstar_column +!use MOM_regridding, only : build_rho_column, build_sigma_column +use MOM_regridding, only : set_regrid_params +use regrid_consts, only : coordinateMode use MOM_verticalGrid, only : verticalGrid_type use diag_axis_mod, only : get_diag_axis_name @@ -104,6 +107,7 @@ module MOM_diag_mediator type, private :: diag_type logical :: in_use integer :: fms_diag_id ! underlying FMS diag id + integer :: vertical_coord type(axes_grp), pointer :: remap_axes => null() real, pointer, dimension(:,:) :: mask2d => null() real, pointer, dimension(:,:,:) :: mask3d => null() @@ -157,23 +161,31 @@ module MOM_diag_mediator ! Needed to diagnostic regridding using ALE type(remapping_CS) :: remap_cs - type(regridding_CS) :: regrid_cs - ! Nominal interface locations for Z remapping + type(regridding_CS) :: zstar_regrid_cs, rho_regrid_cs, sigma_regrid_cs + + ! Nominal interface locations for Z* remapping real, dimension(:), allocatable :: zi_remap - ! Nominal layer locations for Z remapping + ! Nominal layer locations for Z* remapping real, dimension(:), allocatable :: zl_remap - ! Number of z levels used for remapping - integer :: nz_remap + ! Number of vertical z* levels used for remapping + integer :: nz_remap, nrho_remap, nsigma_remap - ! Output grid thicknesses - real, dimension(:,:,:), allocatable :: h_zoutput + ! Remap grid thicknesses for different vertical coords + real, dimension(:,:,:), allocatable :: h_zstar_remap + real, dimension(:,:,:), allocatable :: h_rho_remap + real, dimension(:,:,:), allocatable :: h_sigma_remap ! 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 :: do_z_remap_on_u, do_z_remap_on_v, do_z_remap_on_T + logical :: do_rho_remap_on_u, do_rho_remap_on_v, do_rho_remap_on_T + logical :: do_sigma_remap_on_u, do_sigma_remap_on_v, do_sigma_remap_on_T logical :: remapping_initialized !> String appended to module name for z*-remapped diagnostics character(len=6) :: z_remap_suffix = '_z_new' + !> String appended to module name for rho-remapped diagnostics + character(len=4) :: rho_remap_suffix = '_rho' + character(len=6) :: sigma_remap_suffix = '_sigma' ! Pointer to H and G for remapping real, dimension(:,:,:), pointer :: h => null() @@ -339,7 +351,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) id_zzi = 1 endif - ! Axes for z remapping + ! Axes for vertical coordinate remapping call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL, & x_cell_method='mean', y_cell_method='mean', v_cell_method='mean') call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL, & @@ -691,40 +703,75 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) diag => diag_cs%diags(diag_field_id) do while (associated(diag)) - if (associated(diag%remap_axes)) then - ! Remap this field to another vertical coordinate. + if (.not. associated(diag%remap_axes)) then + ! No remapping + call post_data_3d_low(diag, field, diag_cs, is_static, mask) + else + ! Remap this field to another vertical coordinate. if (present(mask)) then call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") endif - if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) - allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) - call remap_diag_to_z(field, diag, diag_cs, remapped_field) - if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) - if (associated(diag%mask3d)) then - ! Since 3d masks do not vary in the vertical, just use as much as is - ! needed. - call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & - mask=diag%mask3d(:,:,:diag_cs%nz_remap)) - else - call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + if (diag%vertical_coord == coordinateMode('Z*')) then + + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) + allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) + call remap_diag_to_zstar(field, diag, diag_cs, remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) + if (associated(diag%mask3d)) then + ! Since 3d masks do not vary in the vertical, just use as much as is + ! needed. + call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & + mask=diag%mask3d(:,:,:diag_cs%nz_remap)) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) + deallocate(remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) + + elseif (diag%vertical_coord == coordinateMode('RHO')) then + + allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) + call remap_diag_to_rho(field, diag, diag_cs, remapped_field) + if (associated(diag%mask3d)) then + ! Since 3d masks do not vary in the vertical, just use as much as is + ! needed. + call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & + mask=diag%mask3d(:,:,:diag_cs%nz_remap)) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + deallocate(remapped_field) + + elseif (diag%vertical_coord == coordinateMode('SIGMA')) then + + allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) + call remap_diag_to_sigma(field, diag, diag_cs, remapped_field) + if (associated(diag%mask3d)) then + ! Since 3d masks do not vary in the vertical, just use as much as is + ! needed. + call post_data_3d_low(diag, remapped_field, diag_cs, is_static, & + mask=diag%mask3d(:,:,:diag_cs%nz_remap)) + else + call post_data_3d_low(diag, remapped_field, diag_cs, is_static) + endif + deallocate(remapped_field) endif - if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) - deallocate(remapped_field) - if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) - else - call post_data_3d_low(diag, field, diag_cs, is_static, mask) + endif + diag => diag%next enddo + if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) end subroutine post_data_3d !> Remap diagnostic field to z-star vertical grid. !! \note This uses grids generated by diag_update_target_grids() -subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) +subroutine remap_diag_to_zstar(field, diag, diag_cs, remapped_field) real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped type(diag_type), intent(in) :: diag !< Structure containing remapping/masking information type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure @@ -759,7 +806,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo #endif h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) - h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i+1,j,:) ) + h_dest(:) = 0.5 * ( diag_cs%h_zstar_remap(i,j,:) + diag_cs%h_zstar_remap(i+1,j,:) ) call remapping_core_h(nz_src, h_src(:), field(I,j,:), & nz_dest, h_dest(:), remapped_field(I,j,:), diag_cs%remap_cs) do k=1, nz_dest @@ -784,7 +831,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo #endif h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) - h_dest(:) = 0.5 * ( diag_cs%h_zoutput(i,j,:) + diag_cs%h_zoutput(i,j+1,:) ) + h_dest(:) = 0.5 * ( diag_cs%h_zstar_remap(i,j,:) + diag_cs%h_zstar_remap(i,j+1,:) ) call remapping_core_h(nz_src, h_src(:), field(i,J,:), & nz_dest, h_dest(:), remapped_field(i,J,:), diag_cs%remap_cs) do k=1, nz_dest @@ -808,7 +855,7 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) "remap_diag_to_z: H has changed since remapping grids were updated") enddo #endif - h_dest(:) = diag_cs%h_zoutput(i,j,:) + h_dest(:) = diag_cs%h_zstar_remap(i,j,:) call remapping_core_h(nz_src, diag_cs%h(i,j,:), field(i,j,:), & nz_dest, h_dest(:), remapped_field(i,j,:), diag_cs%remap_cs) do k=1, nz_dest @@ -821,7 +868,93 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) enddo endif -end subroutine remap_diag_to_z +end subroutine remap_diag_to_zstar + +!> Remap diagnostic field to isopycnal vertical grid. +!! \note This uses grids generated by diag_update_target_grids() +subroutine remap_diag_to_rho(field, diag, diag_cs, remapped_field) + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped + type(diag_type), intent(in) :: diag !< Structure containing remapping/masking information + type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure + real, dimension(:,:,:), intent(inout) :: remapped_field !< Field argument remapped to new coordinate + ! Local variables + real, dimension(diag_cs%nz_remap+1) :: dz + real, dimension(diag_cs%nz_remap) :: h_dest + real, dimension(size(diag_cs%h, 3)) :: h_src + integer :: nz_src, nz_dest + integer :: i, j, k + + call assert(diag_cs%remapping_initialized, & + 'remap_diag_to_rho: Remmaping not initialized.') + call assert(size(field, 3) == size(diag_cs%h, 3), & + 'remap_diag_to_rho: Remap field and thickness z-axes do not match.') + + remapped_field(:,:,:) = diag_cs%missing_value + nz_src = size(field, 3) + nz_dest = diag_cs%nz_remap + + 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 +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "remap_diag_to_rho: H has changed since remapping grids were updated") + enddo +#endif + h_dest(:) = diag_cs%h_rho_remap(i,j,:) + call remapping_core_h(nz_src, diag_cs%h(i,j,:), field(i,j,:), & + nz_dest, h_dest(:), remapped_field(i,j,:), diag_cs%remap_cs) + enddo + enddo + +end subroutine remap_diag_to_rho + +!> Remap diagnostic field to isopycnal vertical grid. +!! \note This uses grids generated by diag_update_target_grids() +subroutine remap_diag_to_sigma(field, diag, diag_cs, remapped_field) + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped + type(diag_type), intent(in) :: diag !< Structure containing remapping/masking information + type(diag_ctrl), intent(in) :: diag_cs !< Diagnostics control structure + real, dimension(:,:,:), intent(inout) :: remapped_field !< Field argument remapped to new coordinate + ! Local variables + real, dimension(diag_cs%nz_remap+1) :: dz + real, dimension(diag_cs%nz_remap) :: h_dest + real, dimension(size(diag_cs%h, 3)) :: h_src + integer :: nz_src, nz_dest + integer :: i, j, k + + call assert(diag_cs%remapping_initialized, & + 'remap_diag_to_sigma: Remmaping not initialized.') + call assert(size(field, 3) == size(diag_cs%h, 3), & + 'remap_diag_to_sigma: Remap field and thickness z-axes do not match.') + + remapped_field = diag_cs%missing_value + nz_src = size(field, 3) + nz_dest = diag_cs%nz_remap + + 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 +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + if (diag_cs%h_old(i,j,k) /= diag_cs%h(i,j,k)) call MOM_error(FATAL, & + "remap_diag_to_rho: H has changed since remapping grids were updated") + enddo +#endif + h_dest(:) = diag_cs%h_sigma_remap(i,j,:) + call remapping_core_h(nz_src, diag_cs%h(i,j,:), field(i,j,:), & + nz_dest, h_dest(:), remapped_field(i,j,:), diag_cs%remap_cs) + enddo + enddo + +end subroutine remap_diag_to_sigma !> Build/update target vertical grids for diagnostic remapping. !! \note The target grids need to be updated whenever sea surface @@ -835,8 +968,6 @@ subroutine diag_update_target_grids(diag_cs) integer :: nz_src, nz_dest integer :: i, j, k logical :: force, h_changed - real, dimension(size(diag_cs%h, 3)) :: h_src_1d - real, dimension(diag_cs%nz_remap) :: h_zout_1d real, dimension(diag_cs%nz_remap+1) :: z_out_1d nz_dest = diag_cs%nz_remap @@ -855,29 +986,55 @@ subroutine diag_update_target_grids(diag_cs) call assert(allocated(diag_cs%zi_remap), & 'update_diag_target_grids: Remapping axis not initialized') - ! Initialize remapping system, on the first call - call initialize_regridding(nz_dest, 'Z*', 'PPM_IH4', diag_cs%regrid_cs) + ! Initialize remapping and regridding on the first call call initialize_remapping(diag_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false.) - call set_regrid_params(diag_cs%regrid_cs, min_thickness=0.) + + call initialize_regridding(nz_dest, 'Z*', 'PPM_IH4', diag_cs%zstar_regrid_cs) + call set_regrid_params(diag_cs%zstar_regrid_cs, min_thickness=0.) + call setCoordinateResolution(diag_cs%zi_remap(2:) - & + diag_cs%zi_remap(:nz_dest), & + diag_cs%zstar_regrid_cs) + allocate(diag_cs%h_zstar_remap(G%isd:G%ied,G%jsd:G%jed,nz_dest)) + + call initialize_regridding(nz_dest, 'RHO', 'PPM_IH4', diag_cs%rho_regrid_cs) + call set_regrid_params(diag_cs%rho_regrid_cs, min_thickness=0.) + ! FIXME: add set_target_density() + allocate(diag_cs%h_rho_remap(G%isd:G%ied,G%jsd:G%jed,nz_dest)) + + call initialize_regridding(nz_dest, 'SIGMA', 'PPM_IH4', diag_cs%sigma_regrid_cs) + call set_regrid_params(diag_cs%sigma_regrid_cs, min_thickness=0.) + ! FIXME: this needs to be non-dimensional coordinate resolution for sigma, + ! not metres as for zstar. call setCoordinateResolution(diag_cs%zi_remap(2:) - & - diag_cs%zi_remap(:nz_dest), diag_cs%regrid_cs) + diag_cs%zi_remap(:nz_dest), diag_cs%sigma_regrid_cs) + allocate(diag_cs%h_sigma_remap(G%isd:G%ied,G%jsd:G%jed,nz_dest)) - allocate(diag_cs%h_zoutput(G%isd:G%ied,G%jsd:G%jed,nz_dest)) endif - ! Store z-star thicknesses on h-points + ! Calculate remapping thicknesses for different target grids based on + ! nominal/target interface locations. + ! FIXME: are these bounds correct? do j=G%jsc, G%jec+1 do i=G%isc, G%iec+1 if (G%mask2dT(i,j)==0.) then - diag_cs%h_zoutput(i,j,:) = 0. + diag_cs%h_zstar_remap(i,j,:) = 0. + diag_cs%h_rho_remap(i,j,:) = 0. + diag_cs%h_sigma_remap(i,j,:) = 0. cycle endif - h_src_1d(:) = diag_cs%h(i,j,:) - call build_zstar_column(diag_cs%regrid_cs, nz_dest, G%bathyT(i,j), & - sum(h_src_1d(:)), z_out_1d) - h_zout_1d(:) = z_out_1d(1:nz_dest) - z_out_1d(2:nz_dest+1) - diag_cs%h_zoutput(i,j,:) = h_zout_1d(:) + call build_zstar_column(diag_cs%zstar_regrid_cs, nz_dest, G%bathyT(i,j), & + sum(diag_cs%h(i,j,:)), z_out_1d) + diag_cs%h_zstar_remap(i,j,:) = z_out_1d(1:nz_dest) - z_out_1d(2:nz_dest+1) + + !call build_rho_column(diag_cs%rho_regrid_cs, nz_dest, G%bathyT(i,j), & + ! sum(diag_cs%h(i,j,:)), z_out_1d) + !diag_cs%h_rho_remap(i,j,:) = z_out_1d(1:nz_dest) - z_out_1d(2:nz_dest+1) + ! + !call build_sigma_column(diag_cs%zstar_regrid_cs, nz_dest, G%bathyT(i,j), & + ! sum(diag_cs%h(i,j,:)), z_out_1d) + !diag_cs%h_sigma_remap(i,j,:) = z_out_1d(1:nz_dest) - z_out_1d(2:nz_dest+1) + enddo enddo @@ -1075,7 +1232,9 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Local variables real :: MOM_missing_value type(diag_ctrl), pointer :: diag_cs - type(diag_type), pointer :: diag => null(), cmor_diag => null(), z_remap_diag => null(), cmor_z_remap_diag => null() + type(diag_type), pointer :: diag => null(), cmor_diag => null() + type(diag_type), pointer :: z_remap_diag => null(), cmor_z_remap_diag => null() + type(diag_type), pointer :: rho_remap_diag => null(), sigma_remap_diag => null() integer :: primary_id, fms_id, remap_id character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name, cm_string, msg @@ -1087,6 +1246,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & diag => null() cmor_diag => null() z_remap_diag => null() + rho_remap_diag => null() + sigma_remap_diag => null() cmor_z_remap_diag => null() ! Set up the 'primary' diagnostic, first get an underlying FMS id @@ -1112,7 +1273,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & msg, diag_CS, long_name, units, standard_name) endif - ! For the CMOR variation of the both the native and _z diagnostics + ! The CMOR variation if (present(cmor_field_name)) then ! Fallback values for strings set to "NULL" posted_cmor_units = "not provided" ! @@ -1144,7 +1305,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & if (primary_id == -1) then primary_id = get_new_diag_id(diag_cs) endif - ! This will add the cmore variation to the 'primary' diagnostic. + ! This will add the CMOR variation to the 'primary' diagnostic. ! In the case where there is no primary, it will become the primary. call alloc_diag_with_id(primary_id, diag_cs, cmor_diag) cmor_diag%fms_diag_id = fms_id @@ -1171,7 +1332,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) call set_diag_mask(z_remap_diag, diag_cs, axes) - call set_diag_remap_axes(z_remap_diag, diag_cs, axes) + call set_diag_remap_axes(z_remap_diag, diag_cs, axes, 'Z*') + call assert(associated(z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') fms_id = register_diag_field_fms(module_name//trim(diag_cs%z_remap_suffix), field_name, & z_remap_diag%remap_axes%handles, & @@ -1182,13 +1344,13 @@ function register_diag_field(module_name, field_name, axes, init_time, & call attach_cell_methods(fms_id, z_remap_diag%remap_axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method) z_remap_diag%fms_diag_id = fms_id - + if (is_u_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_u = .true. + diag_cs%do_z_remap_on_u = .true. elseif (is_v_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_v = .true. + diag_cs%do_z_remap_on_v = .true. else - diag_cs%do_z_remapping_on_T = .true. + diag_cs%do_z_remap_on_T = .true. endif endif if (is_root_pe() .and. diag_CS%doc_unit > 0) then @@ -1210,7 +1372,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif call alloc_diag_with_id(primary_id, diag_cs, cmor_z_remap_diag) call set_diag_mask(cmor_z_remap_diag, diag_cs, axes) - call set_diag_remap_axes(cmor_z_remap_diag, diag_cs, axes) + call set_diag_remap_axes(cmor_z_remap_diag, diag_cs, axes, 'Z*') call assert(associated(cmor_z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') fms_id = register_diag_field_fms(module_name//trim(diag_cs%z_remap_suffix), cmor_field_name, & cmor_z_remap_diag%remap_axes%handles, & @@ -1221,13 +1383,13 @@ function register_diag_field(module_name, field_name, axes, init_time, & call attach_cell_methods(fms_id, cmor_z_remap_diag%remap_axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method) cmor_z_remap_diag%fms_diag_id = fms_id - + if (is_u_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_u = .true. + diag_cs%do_z_remap_on_u = .true. elseif (is_v_axes(axes, diag_cs)) then - diag_cs%do_z_remapping_on_v = .true. + diag_cs%do_z_remap_on_v = .true. else - diag_cs%do_z_remapping_on_T = .true. + diag_cs%do_z_remap_on_T = .true. endif endif if (is_root_pe() .and. diag_CS%doc_unit > 0) then @@ -1237,6 +1399,50 @@ function register_diag_field(module_name, field_name, axes, init_time, & posted_cmor_standard_name) endif endif + + ! Remap to rho vertical coordinate + if (get_diag_field_id_fms(trim(module_name)//trim(diag_cs%rho_remap_suffix), field_name) /= DIAG_FIELD_NOT_FOUND) then + print*, 'GOT A RHO diagnostic.', trim(module_name)//trim(diag_cs%rho_remap_suffix)//field_name + + if (primary_id == -1) then + primary_id = get_new_diag_id(diag_cs) + endif + call alloc_diag_with_id(primary_id, diag_cs, rho_remap_diag) + call set_diag_mask(rho_remap_diag, diag_cs, axes) + call set_diag_remap_axes(rho_remap_diag, diag_cs, axes, 'RHO') + call assert(associated(rho_remap_diag%remap_axes), 'register_diag_field: remap axes not set') + fms_id = register_diag_field_fms(module_name//trim(diag_cs%rho_remap_suffix), field_name, & + rho_remap_diag%remap_axes%handles, & + init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & + interp_method=interp_method, tile_count=tile_count) + call attach_cell_methods(fms_id, rho_remap_diag%remap_axes, cm_string, & + cell_methods, x_cell_method, y_cell_method, v_cell_method) + rho_remap_diag%fms_diag_id = fms_id + endif + + ! Remap to sigma vertical coordinate + if (get_diag_field_id_fms(trim(module_name)//trim(diag_cs%sigma_remap_suffix), field_name) /= DIAG_FIELD_NOT_FOUND) then + print*, 'GOT A SIGMA diagnostic.', trim(module_name)//trim(diag_cs%sigma_remap_suffix)//field_name + + if (primary_id == -1) then + primary_id = get_new_diag_id(diag_cs) + endif + call alloc_diag_with_id(primary_id, diag_cs, sigma_remap_diag) + call set_diag_mask(sigma_remap_diag, diag_cs, axes) + call set_diag_remap_axes(sigma_remap_diag, diag_cs, axes, 'RHO') + call assert(associated(sigma_remap_diag%remap_axes), 'register_diag_field: remap axes not set') + fms_id = register_diag_field_fms(module_name//trim(diag_cs%sigma_remap_suffix), field_name, & + rho_remap_diag%remap_axes%handles, & + init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & + interp_method=interp_method, tile_count=tile_count) + call attach_cell_methods(fms_id, sigma_remap_diag%remap_axes, cm_string, & + cell_methods, x_cell_method, y_cell_method, v_cell_method) + sigma_remap_diag%fms_diag_id = fms_id + endif endif register_diag_field = primary_id @@ -1684,6 +1890,7 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) do i=1, DIAG_ALLOC_CHUNK_SIZE diag_cs%diags(i)%in_use = .false. diag_cs%diags(i)%next => null() + diag_cs%diags(i)%vertical_coord = coordinateMode('LAYER') diag_cs%diags(i)%remap_axes => null() diag_cs%diags(i)%mask2d => null() diag_cs%diags(i)%mask3d => null() @@ -1693,9 +1900,9 @@ subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) diag_cs%G => G diag_cs%h => null() diag_cs%nz_remap = -1 - 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%do_z_remap_on_u = .false. + diag_cs%do_z_remap_on_v = .false. + diag_cs%do_z_remap_on_T = .false. diag_cs%remapping_initialized = .false. #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,nz)) @@ -1822,7 +2029,7 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) if (allocated(diag_cs%zi_remap)) deallocate(diag_cs%zi_remap) if (allocated(diag_cs%zl_remap)) deallocate(diag_cs%zl_remap) - if (allocated(diag_cs%h_zoutput)) deallocate(diag_cs%h_zoutput) + if (allocated(diag_cs%h_zstar_remap)) deallocate(diag_cs%h_zstar_remap) deallocate(diag_cs%mask3dTL) deallocate(diag_cs%mask3dBuL) deallocate(diag_cs%mask3dCuL) @@ -1862,12 +2069,14 @@ function i2s(a,n_in) i2s = adjustl(i2s) end function i2s -subroutine set_diag_remap_axes(diag, diag_cs, axes) +subroutine set_diag_remap_axes(diag, diag_cs, axes, vertical_coord) ! Associate a remapping axes with the a diagnostic based on the axes info. type(diag_ctrl), target, intent(in) :: diag_cs type(diag_type), pointer, intent(inout) :: diag type(axes_grp), intent(in) :: axes + character(len=*), intent(in) :: vertical_coord + diag%vertical_coord = coordinateMode(vertical_coord) diag%remap_axes => null() if (axes%rank .eq. 3) then if ((axes%id .eq. diag_cs%axesTL%id)) then @@ -2034,6 +2243,8 @@ function get_new_diag_id(diag_cs) ! Initialise new part of the diag array. do i=diag_cs%next_free_diag_id, size(diag_cs%diags) diag_cs%diags(i)%in_use = .false. + diag_cs%diags(i)%vertical_coord = coordinateMode('LAYER') + diag_cs%diags(i)%remap_axes => null() diag_cs%diags(i)%next => null() diag_cs%diags(i)%mask2d => null() diag_cs%diags(i)%mask3d => null()