diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index eac3040df2..087a47ff14 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -37,7 +37,7 @@ module ocean_model_mod use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only : calculate_surface_state use MOM_constants, only : CELSIUS_KELVIN_OFFSET -use MOM_diag_mediator, only : enable_averaging, disable_averaging +use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end use MOM_domains, only : pass_vector, AGRID, BGRID_NE, CGRID_NE use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe @@ -502,6 +502,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) type(ocean_public_type), intent(inout) :: Ocean_sfc type(ocean_state_type), pointer :: Ocean_state type(time_type), intent(in) :: Time + ! This subroutine terminates the model run, saving the ocean state in a ! restart file and deallocating any data associated with the ocean. @@ -512,7 +513,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) ! (in) Time - The model time, used for writing restarts. call ocean_model_save_restart(Ocean_state, Time) - call diag_mediator_end(Time) + call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) end subroutine ocean_model_end diff --git a/config_src/ice_solo_driver/ice_shelf_driver.F90 b/config_src/ice_solo_driver/ice_shelf_driver.F90 index f00071f57b..2e694a2b0a 100644 --- a/config_src/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/ice_solo_driver/ice_shelf_driver.F90 @@ -416,7 +416,7 @@ program SHELF_main close(unit) endif - call diag_mediator_end(Time, end_diag_manager=.true.) + call diag_mediator_end(Time, ice_shelf_CSp%diag, end_diag_manager=.true.) call cpu_clock_end(termClock) call io_infra_end ; call MOM_infra_end diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index f3b72f9415..139adc0255 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -517,7 +517,7 @@ program MOM_main endif call callTree_waypoint("End MOM_main") - call diag_mediator_end(Time, end_diag_manager=.true.) + call diag_mediator_end(Time, MOM_CSp%diag, end_diag_manager=.true.) call cpu_clock_end(termClock) call io_infra_end ; call MOM_infra_end diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 487d12ce8c..704fdad396 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -534,15 +534,22 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) if (k<=n0) then; hTmp = h0(k); else; hTmp = 0.; endif z0 = z0 + hTmp ; z1 = z1 + ( hTmp + ( dx(k+1) - dx(k) ) ) enddo + ! Maximum error based on guess at maximum roundoff if (abs(totalHU2-totalHU0) > (err0+err2)*max(real(n0), real(n1)) .and. (err0+err2)/=0.) then - write(0,*) 'h0=',h0 - write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1) - write(0,*) 'u0=',u0 - write(0,*) 'u1=',u1 - write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0 - write(0,*) 'err0,errF=',err0,err2 - call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//& - 'Total stuff on h0 and hF differ by more than roundoff' ) + ! Maximum relative error + if (abs(totalHU2-totalHU0) / totalHU2 > 1e-09) then + ! Maximum absolute error + if (abs(totalHU2-totalHU0) > 1e-18) then + write(0,*) 'h0=',h0 + write(0,*) 'hf=',h0(1:n1)+dx(2:n1+1)-dx(1:n1) + write(0,*) 'u0=',u0 + write(0,*) 'u1=',u1 + write(0,*) 'total HU0,HUf,f-0=',totalHU0,totalHU2,totalHU2-totalHU0 + write(0,*) 'err0,errF=',err0,err2 + call MOM_error( FATAL, 'MOM_remapping, remapping_core: '//& + 'Total stuff on h0 and hF differ by more than maximum errors' ) + endif + endif endif #endif diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 60b7b9542a..e73b4853ed 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -337,7 +337,8 @@ module MOM use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_coms, only : reproducing_sum -use MOM_diag_mediator, only : diag_mediator_init, enable_averaging, diag_set_thickness_ptr +use MOM_diag_mediator, only : diag_mediator_init, enable_averaging +use MOM_diag_mediator, only : diag_set_thickness_ptr, diag_update_target_grids use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field use MOM_diag_mediator, only : register_scalar_field @@ -977,6 +978,10 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) call cpu_clock_end(id_clock_pass) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) + if (CS%debug) then call uchksum(u,"Post-dia first u", G, haloshift=2) call vchksum(v,"Post-dia first v", G, haloshift=2) @@ -1039,6 +1044,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) else dtth = dt*min(ntstep,n_max-n+1) endif + call enable_averaging(dtth,Time_local+set_time(int(floor(dtth-dt+0.5))), CS%diag) call cpu_clock_begin(id_clock_thick_diff) if (associated(CS%VarMix)) call calc_slope_functions(h, CS%tv, dt, G, CS%VarMix) @@ -1050,6 +1056,11 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call cpu_clock_end(id_clock_pass) call disable_averaging(CS%diag) if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)") + + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) + endif endif @@ -1159,11 +1170,9 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (CS%useMEKE) call step_forward_MEKE(CS%MEKE, h, CS%VarMix%SN_u, CS%VarMix%SN_v, & CS%visc, dt, G, CS%MEKE_CSp) - call disable_averaging(CS%diag) call cpu_clock_end(id_clock_dynamics) - CS%dt_trans = CS%dt_trans + dt if (thermo_does_span_coupling) then do_advection = (CS%dt_trans + 0.5*dt > dt_therm) @@ -1270,7 +1279,12 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call hchksum(CS%tv%S,"Post-ALE S", G, haloshift=1) call check_redundant("Post-ALE ", u, v, G) endif - endif + endif + + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. This needs to + ! happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_uv_T_S_h, G%Domain) @@ -1928,10 +1942,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif call callTree_waypoint("state variables allocated (initialize_MOM)") - ! Set up a pointers h within diag mediator control structure, - ! this needs to occur _after_ CS%h has been allocated. - call diag_set_thickness_ptr(CS%h, diag) - ! Set the fields that are needed for bitwise identical restarting ! the time stepping scheme. call restart_init(G, param_file, CS%restart_CSp) @@ -1980,11 +1990,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") - ! Initialize the diagnostics mask arrays. - ! This step has to be done after call to MOM_initialize_state - ! and before MOM_diagnostics_init - call diag_masks_set(G, CS%missing, diag) - if (CS%use_ALE_algorithm) then ! For now, this has to follow immediately after MOM_initialize_state because ! the call to initialize_ALE can change CS%h, etc. initialize_ALE should @@ -2008,9 +2013,24 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif endif - ! This call sets up the diagnostic axes. - call cpu_clock_begin(id_clock_MOM_init) + ! Initialize the diagnostics mask arrays. + ! This step has to be done after call to MOM_initialize_state + ! and before MOM_diagnostics_init + call diag_masks_set(G, CS%missing, diag) + + ! Set up a pointers h within diag mediator control structure, + ! this needs to occur _after_ CS%h has been allocated. + call diag_set_thickness_ptr(CS%h, diag) + + ! This call sets up the diagnostic axes. These are needed, + ! e.g. to generate the target grids below. call set_axes_info(G, param_file, diag) + + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. This needs to + call diag_update_target_grids(diag) + + call cpu_clock_begin(id_clock_MOM_init) if (CS%use_ALE_algorithm) then call ALE_writeCoordinateFile( CS%ALE_CSp, G, dirs%output_directory ) endif diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 68e4767215..533c8c9248 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -78,7 +78,7 @@ module MOM_dynamics_legacy_split use MOM_diag_mediator, only : diag_mediator_init, enable_averaging use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field -use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids use MOM_domains, only : MOM_domains_init, pass_var, pass_vector use MOM_domains, only : pass_var_start, pass_var_complete use MOM_domains, only : pass_vector_start, pass_vector_complete @@ -631,8 +631,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) call continuity(u, v, h, hp, uh_in, vh_in, dt, G, & - CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, & - visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) + CS%continuity_CSp, OBC=CS%OBC, & + visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, & + BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (BT_cont_BT_thick) then call cpu_clock_begin(id_clock_pass) @@ -716,8 +717,8 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & ! hp = h + dt * div . uh call cpu_clock_begin(id_clock_continuity) call continuity(up, vp, h, hp, uh, vh, dt, G, CS%continuity_CSp, & - CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, CS%visc_rem_v, & - u_av, v_av, BT_cont=CS%BT_cont) + CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, & + CS%visc_rem_v, u_av, v_av, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) @@ -931,6 +932,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av, & uhbt_out, vhbt_out, u, v) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) if (G%nonblocking_updates) then call cpu_clock_begin(id_clock_pass) pid_h = pass_var_start(h, G%Domain) @@ -969,6 +973,9 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call pass_var(h, G%Domain) call cpu_clock_end(id_clock_pass) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index d5d35161ee..1d0f7974f8 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -78,7 +78,7 @@ module MOM_dynamics_split_RK2 use MOM_diag_mediator, only : diag_mediator_init, enable_averaging use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field -use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids use MOM_domains, only : MOM_domains_init use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type @@ -914,6 +914,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & CS%continuity_CSp, CS%uhbt, CS%vhbt, CS%OBC, & CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call do_group_pass(CS%pass_h, G%Domain) call cpu_clock_end(id_clock_pass) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index dd36f3aa19..ff911873f6 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -78,7 +78,7 @@ module MOM_dynamics_unsplit use MOM_diag_mediator, only : diag_mediator_init, enable_averaging use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field -use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl +use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_target_grids use MOM_domains, only : MOM_domains_init, pass_var, pass_vector use MOM_domains, only : pass_var_start, pass_var_complete use MOM_domains, only : pass_vector_start, pass_vector_complete @@ -268,7 +268,8 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & ! uh = u*h ! hp = h + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u, v, h, hp, uh, vh, dt*0.5, G, CS%continuity_CSp, & + OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(hp, G%Domain) @@ -438,6 +439,9 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & call continuity(upp, vpp, hp, h, uh, vh, & (dt*0.5), G, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) call cpu_clock_begin(id_clock_pass) call pass_var(h, G%Domain) call pass_vector(uh, vh, G%Domain) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 57edff0a20..53ca73f4a0 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -282,7 +282,8 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call cpu_clock_begin(id_clock_continuity) ! This is a duplicate caclulation of the last continuity from the previous step ! and could/should be optimized out. -AJA - call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, CS%continuity_CSp, & + OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call cpu_clock_begin(id_clock_pass) call pass_var(hp, G%Domain) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index a6022ab3de..dd0a901e8b 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -38,7 +38,8 @@ module MOM_diag_mediator use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type use MOM_remapping, only : remapping_CS, remapping_core, initialize_remapping, dzFromH1H2 -use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution, buildGridZStarColumn, setRegriddingMinimumThickness +use MOM_regridding, only : regridding_CS, initialize_regridding, setCoordinateResolution +use MOM_regridding, only : buildGridZStarColumn, setRegriddingMinimumThickness use diag_manager_mod, only : diag_manager_init, diag_manager_end use diag_manager_mod, only : send_data, diag_axis_init @@ -58,6 +59,8 @@ module MOM_diag_mediator #define DIM_J(a) lbound(a, 2):ubound(a, 2) #define DIM_K(a) lbound(a, 3):ubound(a, 3) +#define __DO_SAFETY_CHECKS__ + public set_axes_info, post_data, register_diag_field, time_type public post_data_1d_k public safe_alloc_ptr, safe_alloc_alloc @@ -68,6 +71,7 @@ module MOM_diag_mediator public register_scalar_field public defineAxes, diag_masks_set public diag_set_thickness_ptr +public diag_update_target_grids interface post_data module procedure post_data_3d, post_data_2d, post_data_0d @@ -114,7 +118,7 @@ module MOM_diag_mediator type(axesType) :: axesBi, axesTi, axesCui, axesCvi type(axesType) :: axesB1, axesT1, axesCu1, axesCv1 type(axesType) :: axesZi, axesZL - type(axesType) :: axesTzi, axesTZL + type(axesType) :: axesTzi, axesTZL, axesBZL, axesCuZL, axesCvZL ! Mask arrays for diagnostics real, dimension(:,:), pointer :: mask2dT => null() @@ -139,17 +143,35 @@ module MOM_diag_mediator !default missing value to be sent to ALL diagnostics registrations real :: missing_value = 1.0e+20 - ! Interface locations for Z remapping + ! Needed to diagnostic regridding using ALE + type(remapping_CS) :: remap_cs + type(regridding_CS) :: regrid_cs + ! Nominal interface locations for Z remapping real, dimension(:), allocatable :: zi_remap - ! 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 - ! Pointer to H and G for Z remapping + ! 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 remapping real, dimension(:,:,:), pointer :: h => null() type(ocean_grid_type), pointer :: G => null() +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Keep a copy of h so that we know whether it has changed. If it has then + ! need the target grid for vertical remapping needs to have been updated. + real, dimension(:,:,:), allocatable :: h_old +#endif + end type diag_ctrl integer :: doc_unit = -1 @@ -297,8 +319,11 @@ subroutine set_axes_info(G, param_file, diag_cs, set_vertical) id_zzi = 1 endif - ! Axis for z remapping + ! Axes for z remapping call defineAxes(diag_cs, (/ id_xh, id_yh, id_zzl /), diag_cs%axesTZL) + call defineAxes(diag_cs, (/ id_xq, id_yq, id_zzL /), diag_cs%axesBZL) + call defineAxes(diag_cs, (/ id_xq, id_yh, id_zzL /), diag_cs%axesCuZL) + call defineAxes(diag_cs, (/ id_xh, id_yq, id_zzL /), diag_cs%axesCvZL) ! Vertical axes for the interfaces and layers call defineAxes(diag_cs, (/ id_zi /), diag_cs%axesZi) @@ -625,6 +650,8 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) end subroutine post_data_3d subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) +! Remap diagnostic field to z-star vertical grid. +! This uses grids generated by diag_update_target_grids() real, dimension(:,:,:), intent(in) :: field type(diag_type), intent(in) :: diag @@ -636,58 +663,192 @@ subroutine remap_diag_to_z(field, diag, diag_cs, remapped_field) ! (in) diag - structure used to regulate diagnostic output ! (out) remapped_field - the field argument remapped to z star coordinate - type(remapping_CS) :: remap_cs - type(regridding_CS) :: regrid_cs - integer :: nz_src, i, j, k - real, dimension(diag_cs%nz_remap) :: h_new, h_remap - real, dimension(diag_cs%nz_remap+1) :: dz, zi_tmp + 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_z: Remmaping not initialized.') call assert(size(field, 3) == size(diag_cs%h, 3), & - 'Remap field and thickness z-axes do not match.') - call assert(allocated(diag_cs%zi_remap), 'Remapping axis not initialized') + 'remap_diag_to_z: 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 + +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Check that H is up-to-date. + do k=RANGE_K(diag_cs%h) + do j=RANGE_J(diag_cs%h) + do i=RANGE_I(diag_cs%h) + if (diag_cs%h_old(i, j, k) /= diag_cs%h(i, j, k)) then + call MOM_error(FATAL, "remap_diag_to_z: H has changed since "//& + "remapping grids were updated") + endif + enddo + enddo + enddo +#endif - ! Nominal thicknesses to remap to - h_remap(:) = diag_cs%zi_remap(2:) - diag_cs%zi_remap(:diag_cs%nz_remap) - call initialize_regridding(diag_cs%nz_remap, 'Z*', 'PPM_IH4', regrid_cs) - call setRegriddingMinimumThickness(diag_cs%G%Angstrom, regrid_cs) - call setCoordinateResolution(h_remap, regrid_cs) - call initialize_remapping(nz_src, 'PPM_IH4', remap_cs) + if (is_u_axes(diag%remap_axes, diag_cs)) then + 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 - do j=RANGE_J(field) - do i=RANGE_I(field) - if (associated(diag%mask3d)) then - if (diag%mask3d(i,j, 1) == 0.) cycle - endif + h_dest(:) = diag_cs%zi_u(i, j, 2:) - diag_cs%zi_u(i, j, :size(diag_cs%zi_u, 3)-1) + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i+1,j,:)) + call dzFromH1H2(nz_src, h_src(:), nz_dest, h_dest(:), dz) + call remapping_core(diag_cs%remap_cs, nz_src, h_src(:), & + field(i, j, :), nz_dest, dz, remapped_field(i, j, :)) + + ! Lower levels of the remapped data get squashed to follow bathymetry. + ! If their nominal depth is below the bathymetry remove them. + do k=1, nz_dest + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then + remapped_field(i, j, k:nz_dest) = diag_cs%missing_value + exit + endif + enddo + enddo + enddo - ! Calculate thicknesses for new Z* using nominal thicknesses from - ! h_remap, current bathymetry and total thickness. - call buildGridZstarColumn(regrid_cs, diag_cs%nz_remap, diag_cs%G%bathyT(i, j), & - sum(diag_cs%h(i, j, :)), zi_tmp) - ! Calculate how much thicknesses change between source and dest grids, do - ! remapping - zi_tmp = -zi_tmp - h_new(:) = zi_tmp(2:) - zi_tmp(:size(zi_tmp)-1) - call dzFromH1H2(nz_src, diag_cs%h(i, j, :), diag_cs%nz_remap, h_new, dz) - call remapping_core(remap_cs, nz_src, diag_cs%h(i, j, :), field(i,j,:), & - diag_cs%nz_remap, 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, diag_cs%nz_remap - if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then - remapped_field(i, j, k:diag_cs%nz_remap) = diag_cs%missing_value - exit + 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 + + h_dest(:) = diag_cs%zi_v(i, j, 2:) - diag_cs%zi_v(i, j, :size(diag_cs%zi_v, 3)-1) + h_src(:) = 0.5 * (diag_cs%h(i,j,:) + diag_cs%h(i,j+1,:)) + call dzFromH1H2(nz_src, h_src(:), nz_dest, h_dest(:), dz) + call remapping_core(diag_cs%remap_cs, nz_src, h_src(:), & + field(i, j, :), nz_dest, dz, remapped_field(i, j, :)) + do k=1, nz_dest + if (diag_cs%zi_remap(k) >= diag_cs%G%bathyT(i, j)) then + remapped_field(i, j, k:nz_dest) = diag_cs%missing_value + exit + endif + enddo enddo enddo - enddo + else + 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(:) = 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_remap(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 end subroutine remap_diag_to_z +subroutine diag_update_target_grids(diag_cs) +! Build/update target vertical grids for diagnostic remapping. +! +! The target grids need to be updated whenever sea surface +! height changes. + + type(diag_ctrl), intent(inout) :: diag_cs + +! Arguments: +! (inout) diag_cs - structure used to regulate diagnostic output. + + real, dimension(size(diag_cs%h, 3)) :: h_src + type(ocean_grid_type), pointer :: G + real :: depth + integer :: nz_src, nz_dest + integer :: i, j, k + logical :: force, h_changed + + nz_dest = diag_cs%nz_remap + nz_src = size(diag_cs%h, 3) + G => diag_cs%G + + ! The interface positions for z remapping were not provided, so don't do + ! anything. If z remapping of diagnostics is requested then an error will + ! be triggered on post(). See param DIAG_REMAP_Z_GRID_DEF + if (.not. allocated(diag_cs%zi_remap)) then + return + endif + + if (.not. diag_cs%remapping_initialized) then + call assert(allocated(diag_cs%zi_remap), & + 'update_diag_target_grids: Remapping axis not initialized') + + ! Initialise remapping system, on the first call + 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(:nz_dest), diag_cs%regrid_cs) + + allocate(diag_cs%zi_u(G%iscB:G%iecB,G%jsc:G%jec,nz_dest+1)) + allocate(diag_cs%zi_v(G%isc:G%iec,G%jscB:G%jecB,nz_dest+1)) + allocate(diag_cs%zi_T(G%isc:G%iec,G%jsc:G%jec,nz_dest+1)) + endif + + ! Build z-star grid on u points + if (diag_cs%do_z_remapping_on_u .or. (.not. diag_cs%remapping_initialized)) then + do j=G%jsc, G%jec + do i=G%iscB, G%iecB + 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(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 .or. (.not. diag_cs%remapping_initialized)) then + 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(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 .or. (.not. diag_cs%remapping_initialized)) then + do j=G%jsc, G%jec + do i=G%isc, G%iec + 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 + + diag_cs%remapping_initialized = .true. +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + ! Keep a copy of H - used to check whether grids are up-to-date + ! when doing remapping. + diag_cs%h_old(:,:,:) = diag_cs%h(:,:,:) +#endif + +end subroutine diag_update_target_grids + subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) type(diag_type), intent(in) :: diag real, intent(in) :: field(:,:,:) @@ -946,30 +1107,39 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif endif - ! Remap to z vertical coordinate, note that only - ! diagnostics on T points and layers (not interfaces) are supported, - ! for other diagnostics the old remapping approach can still be used - if (axes%id == diag_cs%axesTL%id) then - fms_id = register_diag_field_fms(module_name//'_z_new', field_name, diag_cs%axesTZL%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) - if (fms_id /= DIAG_FIELD_NOT_FOUND) then - ! Check that we have read in necessary remapping information from file - if (.not. allocated(diag_cs%zi_remap)) then - call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & - 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') - endif + ! Remap to z vertical coordinate, note that only diagnostics on layers + ! (not interfaces) are supported, also B axes are not supported yet + if (get_diag_field_id_fms(module_name//'_z_new', field_name) /= DIAG_FIELD_NOT_FOUND & + .and. is_layer_axes(axes, diag_cs) .and. (.not. is_B_axes(axes, diag_cs)) & + .and. axes%rank == 3) then + if (.not. allocated(diag_cs%zi_remap)) then + call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & + 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') + endif + if (primary_id == -1) then + primary_id = get_new_diag_id(diag_cs) + 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 assert(associated(z_remap_diag%remap_axes), & + 'register_diag_field: remap axes not set') + fms_id = register_diag_field_fms(module_name//'_z_new', field_name, & + z_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) + z_remap_diag%fms_diag_id = fms_id - if (primary_id == -1) then - primary_id = get_new_diag_id(diag_cs) - endif - call alloc_diag_with_id(primary_id, diag_cs, z_remap_diag) - z_remap_diag%fms_diag_id = fms_id - z_remap_diag%remap_axes => diag_cs%axesTZL - call set_diag_mask(z_remap_diag, diag_cs, diag_cs%axesTZL) + if (is_u_axes(axes, diag_cs)) then + diag_cs%do_z_remapping_on_u = .true. + elseif (is_v_axes(axes, diag_cs)) then + diag_cs%do_z_remapping_on_v = .true. + else + diag_cs%do_z_remapping_on_T = .true. endif + endif ! Document diagnostics in list of available diagnostics @@ -981,7 +1151,8 @@ function register_diag_field(module_name, field_name, axes, init_time, & posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif - if (axes%id == diag_cs%axesTL%id) then + if (is_layer_axes(axes, diag_cs) .and. (.not. is_B_axes(axes, diag_cs)) & + .and. axes%rank == 3) then call log_available_diag(associated(z_remap_diag), module_name//'_z_new', field_name, & long_name, units, standard_name) endif @@ -1327,7 +1498,14 @@ 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%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%remapping_initialized = .false. + allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,G%ks:G%ke)) + diag_cs%h_old(:,:,:) = 0.0 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) @@ -1426,14 +1604,36 @@ subroutine diag_mediator_close_registration() end subroutine diag_mediator_close_registration -subroutine diag_mediator_end(time, end_diag_manager) - type(time_type), intent(in) :: time - logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end() +subroutine diag_mediator_end(time, diag_cs, end_diag_manager) + type(time_type), intent(in) :: time + type(diag_ctrl), intent(inout) :: diag_cs + logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end() if (doc_unit > -1) then close(doc_unit) ; doc_unit = -3 endif + deallocate(diag_cs%diags) + + 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%zi_u)) deallocate(diag_cs%zi_u) + if (allocated(diag_cs%zi_v)) deallocate(diag_cs%zi_v) + if (allocated(diag_cs%zi_T)) deallocate(diag_cs%zi_T) + deallocate(diag_cs%mask3dTL) + deallocate(diag_cs%mask3dBuL) + deallocate(diag_cs%mask3dCuL) + deallocate(diag_cs%mask3dCvL) + deallocate(diag_cs%mask3dTi) + deallocate(diag_cs%mask3dBui) + deallocate(diag_cs%mask3dCui) + deallocate(diag_cs%mask3dCvi) + +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + deallocate(diag_cs%h_old) +#endif + if (present(end_diag_manager)) then if (end_diag_manager) call diag_manager_end(time) endif @@ -1460,18 +1660,39 @@ function i2s(a,n_in) i2s = adjustl(i2s) end function i2s +subroutine set_diag_remap_axes(diag, diag_cs, axes) + ! Associate a remapping axes with the a diagnostic based on the axes info. + type(diag_ctrl), target, intent(inout) :: diag_cs + type(diag_type), pointer, intent(out) :: diag + type(axesType), intent(in) :: axes + + diag%remap_axes => null() + if (axes%rank .eq. 3) then + if ((axes%id .eq. diag_cs%axesTL%id)) then + diag%remap_axes => diag_cs%axesTZL + elseif(axes%id .eq. diag_cs%axesBL%id) then + diag%remap_axes => diag_cs%axesBZL + elseif(axes%id .eq. diag_cs%axesCuL%id ) then + diag%remap_axes => diag_cs%axesCuZL + elseif(axes%id .eq. diag_cs%axesCvL%id) then + diag%remap_axes => diag_cs%axesCvZL + endif + endif + +end subroutine set_diag_remap_axes + subroutine set_diag_mask(diag, diag_cs, axes) ! Associate a mask with the a diagnostic based on the axes info. - type(diag_ctrl), intent(inout) :: diag_cs + type(diag_ctrl), target, intent(inout) :: diag_cs type(diag_type), pointer, intent(out) :: diag type(axesType), intent(in) :: axes + diag%mask2d => null() + diag%mask3d => null() + if (axes%rank .eq. 3) then - diag%mask2d => null() - diag%mask3d => null() - if ((axes%id .eq. diag_cs%axesTL%id) .or. & - (axes%id .eq. diag_cs%axesTZL%id)) then + if ((axes%id .eq. diag_cs%axesTL%id)) then diag%mask3d => diag_cs%mask3dTL elseif(axes%id .eq. diag_cs%axesBL%id) then diag%mask3d => diag_cs%mask3dBuL @@ -1491,8 +1712,6 @@ subroutine set_diag_mask(diag, diag_cs, axes) !call assert(associated(diag%mask3d), "MOM_diag_mediator.F90: Invalid 3d axes id") elseif(axes%rank .eq. 2) then - diag%mask2d => null() - diag%mask3d => null() if (axes%id .eq. diag_cs%axesT1%id) then diag%mask2d => diag_cs%mask2dT elseif(axes%id .eq. diag_cs%axesB1%id) then @@ -1508,6 +1727,83 @@ subroutine set_diag_mask(diag, diag_cs, axes) end subroutine set_diag_mask +function is_layer_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_layer_axes + + is_layer_axes = .false. + + if (axes%id == diag_cs%axesTZL%id .or. & + axes%id == diag_cs%axesBZL%id .or. & + axes%id == diag_cs%axesCuZL%id .or. & + axes%id == diag_cs%axesCvZL%id .or. & + axes%id == diag_cs%axesZL%id .or. & + axes%id == diag_cs%axesTL%id .or. & + axes%id == diag_cs%axesBL%id .or. & + axes%id == diag_cs%axesCuL%id .or. & + axes%id == diag_cs%axesCvL%id) then + is_layer_axes = .true. + endif + +end function is_layer_axes + +function is_u_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_u_axes + + is_u_axes = .false. + + if (axes%id == diag_cs%axesCuZL%id .or. & + axes%id == diag_cs%axesCuL%id .or. & + axes%id == diag_cs%axesCui%id .or. & + axes%id == diag_cs%axesCu1%id) then + is_u_axes = .true. + endif + +end function is_u_axes + +function is_v_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_v_axes + + is_v_axes = .false. + + if (axes%id == diag_cs%axesCuZL%id .or. & + axes%id == diag_cs%axesCuL%id .or. & + axes%id == diag_cs%axesCui%id .or. & + axes%id == diag_cs%axesCu1%id) then + is_v_axes = .true. + endif + +end function is_v_axes + +function is_B_axes(axes, diag_cs) + + type(axesType), intent(in) :: axes + type(diag_ctrl), intent(in) :: diag_cs + + logical :: is_B_axes + + is_B_axes = .false. + + if (axes%id == diag_cs%axesBZL%id .or. & + axes%id == diag_cs%axesBL%id .or. & + axes%id == diag_cs%axesBi%id .or. & + axes%id == diag_cs%axesB1%id) then + is_B_axes = .true. + endif + +end function is_B_axes + ! Allocate a new diagnostic id, it may be necessary to expand the diagnostics ! array. function get_new_diag_id(diag_cs) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index bae71a5933..41286a7a87 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -61,6 +61,7 @@ module MOM_mixed_layer_restrat use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type +use MOM_diag_mediator, only : diag_update_target_grids use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_forcing_type, only : forcing @@ -394,6 +395,11 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, G, CS) enddo ; enddo ; enddo !$OMP end parallel + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + ! This needs to happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) + ! Offer fields for averaging. if (query_averaging_enabled(CS%diag)) then if (CS%id_urestrat_time > 0) call post_data(CS%id_urestrat_time, utimescale_diag, CS%diag) @@ -644,6 +650,10 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, fluxes, dt, G, CS) enddo ; enddo ; enddo !$OMP end parallel + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) + ! Offer fields for averaging. if (query_averaging_enabled(CS%diag) .and. & ((CS%id_urestrat_time > 0) .or. (CS%id_vrestrat_time > 0))) then diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index ab936f2a5e..91d95bee15 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -48,6 +48,7 @@ module MOM_thickness_diffuse use MOM_checksums, only : hchksum, uchksum, vchksum use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type +use MOM_diag_mediator, only : diag_update_target_grids use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_EOS, only : calculate_density, calculate_density_derivs use MOM_file_parser, only : get_param, log_version, param_file_type @@ -339,6 +340,11 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CDp, CS) enddo ; enddo enddo + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + ! This needs to happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) + if (MEKE_not_null .AND. ASSOCIATED(VarMix)) then if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then !$OMP parallel do default(none) shared(is,ie,js,je,MEKE,VarMix) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index da224979ca..a1313727b3 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -55,7 +55,7 @@ module MOM_bulk_mixed_layer use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc -use MOM_diag_mediator, only : time_type, diag_ctrl +use MOM_diag_mediator, only : time_type, diag_ctrl, diag_update_target_grids use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -805,6 +805,11 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, CS, & enddo ! j loop + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + ! This needs to happen after the H update and before the next post_data. + call diag_update_target_grids(CS%diag) + !$OMP end parallel if (write_diags) then diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 1f64ab5e9a..b408a0635b 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -72,7 +72,7 @@ module MOM_diabatic_driver use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr -use MOM_diag_mediator, only : diag_ctrl, time_type +use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_target_grids use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags use MOM_diffConvection, only : diffConvection_CS, diffConvection_init use MOM_diffConvection, only : diffConvection_calculate, diffConvection_end @@ -410,6 +410,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G) endif + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) + ! Set_opacity estimates the optical properties of the water column. ! It will need to be modified later to include information about the ! biological properties and layer thicknesses. @@ -444,6 +448,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) G, CS%bulkmixedlayer_CSp, CS%optics, & CS%aggregate_FW_forcing, dt, last_call=.true.) endif + ! Keep salinity from falling below a small but positive threshold. ! This occurs when the ice model attempts to extract more salt than ! is actually present in the ocean. @@ -912,6 +917,10 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) endif + ! Whenever thickness changes let the diag manager know, target grids + ! for vertical remapping may need to be regenerated. + call diag_update_target_grids(CS%diag) + if ((CS%id_Tdif > 0) .or. (CS%id_Tdif_z > 0) .or. & (CS%id_Tadv > 0) .or. (CS%id_Tadv_z > 0)) then do j=js,je ; do i=is,ie