diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index ca3b9d54de..8116ba3e17 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -10,9 +10,9 @@ module MOM_ALE ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_debugging, only : check_column_integrals, hchksum, uvchksum +use MOM_debugging, only : check_column_integrals use MOM_diag_mediator, only : register_diag_field, post_data, diag_ctrl -use MOM_diag_mediator, only : time_type, diag_update_remap_grids +use MOM_diag_mediator, only : time_type, diag_update_remap_grids, query_averaging_enabled use MOM_diag_vkernels, only : interpolate_column, reintegrate_column use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, FATAL, WARNING @@ -123,12 +123,12 @@ module MOM_ALE ! Publicly available functions public ALE_init public ALE_end -public ALE_main -public ALE_main_offline +public ALE_regrid public ALE_offline_inputs -public ALE_offline_tracer_final public ALE_regrid_accelerated public ALE_remap_scalar +public ALE_remap_tracers +public ALE_remap_velocities public ALE_PLM_edge_values public TS_PLM_edge_values public TS_PPM_edge_values @@ -140,6 +140,8 @@ module MOM_ALE public ALE_updateVerticalGridType public ALE_initThicknessToCoord public ALE_update_regrid_weights +public pre_ALE_diagnostics +public pre_ALE_adjustments public ALE_remap_init_conds public ALE_register_diags @@ -164,7 +166,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) ! Local variables character(len=40) :: mdl = "MOM_ALE" ! This module's name. character(len=80) :: string, vel_string ! Temporary strings - real :: filter_shallow_depth, filter_deep_depth + real :: filter_shallow_depth, filter_deep_depth ! Depth ranges of filtering [H ~> m or kg m-2] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. logical :: answers_2018 ! If true, use the order of arithmetic and expressions for remapping @@ -392,11 +394,9 @@ subroutine ALE_end(CS) end subroutine ALE_end -!> Takes care of (1) building a new grid and (2) remapping all variables between -!! the old grid and the new grid. The creation of the new grid can be based -!! on z coordinates, target interface densities, sigma coordinates or any -!! arbitrary coordinate system. -subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) +!> Save any diagnostics of the state before ALE remapping. These diagnostics are +!! mostly used for debugging. +subroutine pre_ALE_diagnostics(G, GV, US, h, u, v, tv, CS) type(ocean_grid_type), intent(in) :: G !< Ocean grid informations type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -405,23 +405,11 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocity field [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< Meridional velocity field [L T-1 ~> m s-1] type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure - type(tracer_registry_type), pointer :: Reg !< Tracer registry structure type(ALE_CS), pointer :: CS !< Regridding parameters and options - type(ocean_OBC_type), pointer :: OBC !< Open boundary structure - real, optional, intent(in) :: dt !< Time step between calls to ALE_main [T ~> s] - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim] - ! Local variables - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg m-2] - logical :: PCM_cell(SZI_(G),SZJ_(G),SZK_(GV)) !< If true, PCM remapping should be used in a cell. - integer :: ntr, i, j, k, isc, iec, jsc, jec, nk - - isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke - if (CS%show_call_tree) call callTree_enter("ALE_main(), MOM_ALE.F90") + ! Local variables + real :: eta_preale(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights before remapping [Z ~> m] - ! These diagnostics of the state before ALE is applied are mostly used for debugging. if (CS%id_u_preale > 0) call post_data(CS%id_u_preale, u, CS%diag) if (CS%id_v_preale > 0) call post_data(CS%id_v_preale, v, CS%diag) if (CS%id_h_preale > 0) call post_data(CS%id_h_preale, h, CS%diag) @@ -432,121 +420,77 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) call post_data(CS%id_e_preale, eta_preale, CS%diag) endif - if (present(dt)) then - call ALE_update_regrid_weights( dt, CS ) - endif - dzRegrid(:,:,:) = 0.0 - - ! If necessary, do some preparatory work to clean up the model state before regridding. +end subroutine pre_ALE_diagnostics - ! This adjusts the input thicknesses prior to remapping, based on the verical coordinate. - if (CS%do_conv_adj) call convective_adjustment(G, GV, h, tv) - if (CS%use_hybgen_unmix) then - ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr - call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h) - endif - ! Build new grid. The new grid is stored in h_new. The old grid is h. - ! Both are needed for the subsequent remapping of variables. - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, conv_adjust=.false., & - frac_shelf_h=frac_shelf_h, PCM_cell=PCM_cell) - - call check_grid( G, GV, h, 0. ) - - if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_main)") - - ! The presence of dt is used for expediency to distinguish whether ALE_main is being called during init - ! or in the main loop. Tendency diagnostics in remap_all_state_vars also rely on this logic. - if (present(dt)) then - call diag_update_remap_grids(CS%diag) - endif - - ! Remap all variables from old grid h onto new grid h_new - call remap_all_state_vars( CS, G, GV, h, h_new, Reg, OBC, dzRegrid, u, v, & - CS%show_call_tree, dt, PCM_cell=PCM_cell ) - - if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_main)") - - ! Override old grid with new one. The new grid 'h_new' is built in - ! one of the 'build_...' routines above. - !$OMP parallel do default(shared) - do k=1,nk ; do j=jsc-1,jec+1 ; do i=isc-1,iec+1 - h(i,j,k) = h_new(i,j,k) - enddo ; enddo ; enddo - - if (CS%debug) then - call hchksum(h, "Post-ALE_main h", G%HI, haloshift=0, scale=GV%H_to_m) - call hchksum(tv%T, "Post-ALE_main T", G%HI, haloshift=0, scale=US%C_to_degC) - call hchksum(tv%S, "Post-ALE_main S", G%HI, haloshift=0, scale=US%S_to_ppt) - call uvchksum("Post-ALE_main [uv]", u, v, G%HI, haloshift=0, scale=US%L_T_to_m_s) - endif - - if (CS%id_dzRegrid>0 .and. present(dt)) call post_data(CS%id_dzRegrid, dzRegrid, CS%diag) - - if (CS%show_call_tree) call callTree_leave("ALE_main()") - -end subroutine ALE_main - -!> Takes care of (1) building a new grid and (2) remapping all variables between -!! the old grid and the new grid. The creation of the new grid can be based -!! on z coordinates, target interface densities, sigma coordinates or any -!! arbitrary coordinate system. -subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, OBC, dt) +!> Potentially do some preparatory work, such as convective adjustment, to clean up the model +!! state before regridding. +subroutine pre_ALE_adjustments(G, GV, US, h, tv, Reg, CS, u, v) type(ocean_grid_type), intent(in) :: G !< Ocean grid informations type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the !! last time step [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure type(tracer_registry_type), pointer :: Reg !< Tracer registry structure type(ALE_CS), pointer :: CS !< Regridding parameters and options - type(ocean_OBC_type), pointer :: OBC !< Open boundary structure - real, optional, intent(in) :: dt !< Time step between calls to ALE_main [T ~> s] - ! Local variables - real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg m-2] - integer :: ntr, i, j, k, isc, iec, jsc, jec, nk - - isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke - - if (CS%show_call_tree) call callTree_enter("ALE_main_offline(), MOM_ALE.F90") + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + optional, intent(inout) :: u !< Zonal velocity field [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + optional, intent(inout) :: v !< Meridional velocity field [L T-1 ~> m s-1] - if (present(dt)) then - call ALE_update_regrid_weights( dt, CS ) - endif - dzRegrid(:,:,:) = 0.0 + integer :: ntr - ! This adjusts the input state prior to remapping, depending on the verical coordinate. + ! Do column-wise convective adjustment. + ! Tracers and velocities should probably also undergo consistent adjustments. if (CS%do_conv_adj) call convective_adjustment(G, GV, h, tv) + if (CS%use_hybgen_unmix) then ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr - call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h) + call hybgen_unmix(G, GV, US, CS%hybgen_unmixCS, tv, Reg, ntr, h) endif - ! Build new grid. The new grid is stored in h_new. The old grid is h. - ! Both are needed for the subsequent remapping of variables. - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, conv_adjust=.false. ) +end subroutine pre_ALE_adjustments - call check_grid( G, GV, h, 0. ) +!> Takes care of building a new grid. The creation of the new grid can be based on z coordinates, +!! target interface densities, sigma coordinates or any arbitrary coordinate system. +subroutine ALE_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_cell) + type(ocean_grid_type), intent(in) :: G !< Ocean grid informations + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses in 3D grid before + !! regridding [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_new !< Layer thicknesses in 3D grid after + !! regridding [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: dzRegrid !< The change in grid interface positions + !! due to regridding, in the same units as + !! thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure + type(ALE_CS), pointer :: CS !< Regridding parameters and options + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim] + logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + optional, intent(out) :: PCM_cell !< If true, use PCM remapping in a cell. - if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_main)") + ! Local variables + logical :: showCallTree - ! Remap all variables from old grid h onto new grid h_new + showCallTree = callTree_showQuery() - call remap_all_state_vars( CS, G, GV, h, h_new, Reg, OBC, debug=CS%show_call_tree, dt=dt ) + if (showCallTree) call callTree_enter("ALE_regrid(), MOM_ALE.F90") - if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_main)") + ! Build the new grid and store it in h_new. The old grid is retained as h. + ! Both are needed for the subsequent remapping of variables. + dzRegrid(:,:,:) = 0.0 + call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, & + frac_shelf_h=frac_shelf_h, PCM_cell=PCM_cell) - ! Override old grid with new one. The new grid 'h_new' is built in - ! one of the 'build_...' routines above. - !$OMP parallel do default(shared) - do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1 - h(i,j,k) = h_new(i,j,k) - enddo ; enddo ; enddo + if (CS%id_dzRegrid>0) then ; if (query_averaging_enabled(CS%diag)) then + call post_data(CS%id_dzRegrid, dzRegrid, CS%diag, alt_h=h_new) + endif ; endif - if (CS%show_call_tree) call callTree_leave("ALE_main()") - if (CS%id_dzRegrid>0 .and. present(dt)) call post_data(CS%id_dzRegrid, dzRegrid, CS%diag) + if (showCallTree) call callTree_leave("ALE_regrid()") -end subroutine ALE_main_offline +end subroutine ALE_regrid !> Regrid/remap stored fields used for offline tracer integrations. These input fields are assumed to have !! the same layer thicknesses at the end of the last offline interval (which should be a Zstar grid). This @@ -560,15 +504,15 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) type(tracer_registry_type), pointer :: Reg !< Tracer registry structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Zonal mass fluxes [H L2 ~> m3 or kg] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Meridional mass fluxes [H L2 ~> m3 or kg] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kd !< Input diffusivites [Z2 T-1 ~> m2 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kd !< Input diffusivities [Z2 T-1 ~> m2 s-1] logical, intent(in ) :: debug !< If true, then turn checksums type(ocean_OBC_type), pointer :: OBC !< Open boundary structure ! Local variables integer :: nk, i, j, k, isc, iec, jsc, jec real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! Layer thicknesses after regridding [H ~> m or kg m-2] - real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: h_src ! Source grid thicknesses at velocity points [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: h_dest ! Destination grid thicknesses at velocity points [H ~> m or kg m-2] + real, dimension(SZK_(GV)) :: h_dest ! Destination grid thicknesses at velocity points [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: temp_vec ! Transports on the destination grid [H L2 ~> m3 or kg] isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke @@ -580,12 +524,11 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) ! Build new grid from the Zstar state onto the requested vertical coordinate. The new grid is stored ! in h_new. The old grid is h. Both are needed for the subsequent remapping of variables. Convective ! adjustment right now is not used because it is unclear what to do with vanished layers - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, conv_adjust=.false. ) - call check_grid( G, GV, h_new, 0. ) + call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid) if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_offline_inputs)") ! Remap all variables from old grid h onto new grid h_new - call remap_all_state_vars( CS, G, GV, h, h_new, Reg, OBC, debug=CS%show_call_tree ) + call ALE_remap_tracers(CS, G, GV, h, h_new, Reg, debug=CS%show_call_tree) if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_inputs)") ! Reintegrate mass transports from Zstar to the offline vertical coordinate @@ -622,93 +565,13 @@ subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC) ! Copy over the new layer thicknesses do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1 - h(i,j,k) = h_new(i,j,k) + h(i,j,k) = h_new(i,j,k) enddo ; enddo ; enddo if (CS%show_call_tree) call callTree_leave("ALE_offline_inputs()") end subroutine ALE_offline_inputs -!> Remaps all tracers from h onto h_target. This is intended to be called when tracers -!! are done offline. In the case where transports don't quite conserve, we still want to -!! make sure that layer thicknesses offline do not drift too far away from the online model -subroutine ALE_offline_tracer_final( G, GV, h, tv, h_target, Reg, CS, OBC) - type(ocean_grid_type), intent(in) :: G !< Ocean grid informations - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the - !! last time step [H ~> m or kg m-2] - type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_target !< Current 3D grid obtained after - !! last time step [H ~> m or kg m-2] - type(tracer_registry_type), pointer :: Reg !< Tracer registry structure - type(ALE_CS), pointer :: CS !< Regridding parameters and options - type(ocean_OBC_type), pointer :: OBC !< Open boundary structure - ! Local variables - - real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid !< The change in grid interface positions - real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new !< Regridded target thicknesses - integer :: ntr, i, j, k, isc, iec, jsc, jec, nk - - isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke - - if (CS%show_call_tree) call callTree_enter("ALE_offline_tracer_final(), MOM_ALE.F90") - ! Need to make sure that h_target is consistent with the current offline ALE confiuration - if (CS%do_conv_adj) call convective_adjustment(G, GV, h_target, tv) - if (CS%use_hybgen_unmix) then - ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr - call hybgen_unmix(G, GV, G%US, CS%hybgen_unmixCS, tv, Reg, ntr, h) - endif - call regridding_main( CS%remapCS, CS%regridCS, G, GV, h_target, tv, h_new, dzRegrid, conv_adjust=.false. ) - call check_grid( G, GV, h_target, 0. ) - - - if (CS%show_call_tree) call callTree_waypoint("Source and target grids checked (ALE_offline_tracer_final)") - - ! Remap all variables from old grid h onto new grid h_new - - call remap_all_state_vars( CS, G, GV, h, h_new, Reg, OBC, debug=CS%show_call_tree ) - - if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_offline_tracer_final)") - - ! Override old grid with new one. The new grid 'h_new' is built in - ! one of the 'build_...' routines above. - !$OMP parallel do default(shared) - do k = 1,nk - do j = jsc-1,jec+1 ; do i = isc-1,iec+1 - h(i,j,k) = h_new(i,j,k) - enddo ; enddo - enddo - if (CS%show_call_tree) call callTree_leave("ALE_offline_tracer_final()") -end subroutine ALE_offline_tracer_final - -!> Check grid for negative thicknesses -subroutine check_grid( G, GV, h, threshold ) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Current 3D grid obtained after the - !! last time step [H ~> m or kg m-2] - real, intent(in) :: threshold !< Value below which to flag issues, - !! [H ~> m or kg m-2] - ! Local variables - integer :: i, j - - do j = G%jsc,G%jec ; do i = G%isc,G%iec - if (G%mask2dT(i,j)>0.) then - if (minval(h(i,j,:)) < threshold) then - write(0,*) 'check_grid: i,j=',i,j,'h(i,j,:)=',h(i,j,:) - if (threshold <= 0.) then - call MOM_error(FATAL,"MOM_ALE, check_grid: negative thickness encountered.") - else - call MOM_error(FATAL,"MOM_ALE, check_grid: too tiny thickness encountered.") - endif - endif - endif - enddo ; enddo - - -end subroutine check_grid - - !> For a state-based coordinate, accelerate the process of regridding by !! repeatedly applying the grid calculation algorithm subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, dzRegrid, initial) @@ -728,7 +591,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d optional, pointer :: Reg !< Tracer registry to remap onto new grid real, optional, intent(in) :: dt !< Model timestep to provide a timescale for regridding [T ~> s] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & - optional, intent(inout) :: dzRegrid !< Final change in interface positions + optional, intent(inout) :: dzRegrid !< Final change in interface positions [H ~> m or kg m-2] logical, optional, intent(in) :: initial !< Whether we're being called from an initialization !! routine (and expect diagnostics to work) @@ -736,11 +599,15 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d integer :: i, j, itt, nz type(thermo_var_ptrs) :: tv_local ! local/intermediate temp/salt type(group_pass_type) :: pass_T_S_h ! group pass if the coordinate has a stencil - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc, h_orig ! A working copy of layer thicknesses - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: T, S ! local temporary state + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc ! A working copy of layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_orig ! The original layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: T ! local temporary temperatures [C ~> degC] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: S ! local temporary salinities [S ~> ppt] ! we have to keep track of the total dzInterface if for some reason ! we're using the old remapping algorithm for u/v - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzInterface, dzIntTotal + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzInterface ! Interface height changes within + ! an iteration [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzIntTotal ! Cumulative interface position changes [H ~> m or kg m-2] real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] nz = GV%ke @@ -775,14 +642,14 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 endif - do itt = 1, n_itt call do_group_pass(pass_T_S_h, G%domain) ! generate new grid if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local) - call regridding_main(CS%remapCS, CS%regridCS, G, GV, h_loc, tv_local, h, dzInterface, conv_adjust=.false.) + + call regridding_main(CS%remapCS, CS%regridCS, G, GV, h_loc, tv_local, h, dzInterface) dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:) ! remap from original grid onto new grid @@ -798,20 +665,18 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d enddo ! remap all state variables (including those that weren't needed for regridding) - call remap_all_state_vars(CS, G, GV, h_orig, h, Reg, OBC, dzIntTotal, u, v) + call ALE_remap_tracers(CS, G, GV, h_orig, h, Reg) + call ALE_remap_velocities(CS, G, GV, h_orig, h, u, v, OBC, dzIntTotal) ! save total dzregrid for diags if needed? if (present(dzRegrid)) dzRegrid(:,:,:) = dzIntTotal(:,:,:) end subroutine ALE_regrid_accelerated -!> This routine takes care of remapping all variable between the old and the -!! new grids. When velocity components need to be remapped, thicknesses at -!! velocity points are taken to be arithmetic averages of tracer thicknesses. -!! This routine is called during initialization of the model at time=0, to +!> This routine takes care of remapping all tracer variables between the old and the +!! new grids. This routine is called during initialization of the model at time=0, to !! remap initial conditions to the model grid. It is also called during a !! time step to update the state. -subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & - dzInterface, u, v, debug, dt, PCM_cell) +subroutine ALE_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell) type(ALE_CS), intent(in) :: CS !< ALE control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -820,25 +685,13 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid !! [H ~> m or kg m-2] type(tracer_registry_type), pointer :: Reg !< Tracer registry structure - type(ocean_OBC_type), pointer :: OBC !< Open boundary structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & - optional, intent(in) :: dzInterface !< Change in interface position - !! [H ~> m or kg m-2] - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - optional, intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - optional, intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1] logical, optional, intent(in) :: debug !< If true, show the call tree real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s] logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: PCM_cell !< Use PCM remapping in cells where true ! Local variables - real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! The vertically summed thicknesses [H ~> m or kg m-2] - real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2] - real, dimension(GV%ke+1) :: dz ! The change in interface heights interpolated to - ! a velocity point [H ~> m or kg m-2] - real :: tr_column(GV%ke) ! A column of updated tracer concentrations + real :: tr_column(GV%ke) ! A column of updated tracer concentrations [CU ~> Conc] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_conc ! The rate of change of concentrations [Conc T-1 ~> Conc s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_cont ! The rate of change of cell-integrated tracer ! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1] or @@ -847,10 +700,6 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & ! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1] logical :: PCM(GV%ke) ! If true, do PCM remapping from a cell. real :: Idt ! The inverse of the timestep [T-1 ~> s-1] - real :: u_src(GV%ke) ! A column of u-velocities on the source grid [L T-1 ~> m s-1] - real :: u_tgt(GV%ke) ! A column of u-velocities on the target grid [L T-1 ~> m s-1] - real :: v_src(GV%ke) ! A column of v-velocities on the source grid [L T-1 ~> m s-1] - real :: v_tgt(GV%ke) ! A column of v-velocities on the target grid [L T-1 ~> m s-1] real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] @@ -861,13 +710,6 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & show_call_tree = .false. if (present(debug)) show_call_tree = debug - ! If remap_uv_using_old_alg is .true. and u or v is requested, then we must have dzInterface. Otherwise, - ! u and v can be remapped without dzInterface - if ( .not. present(dzInterface) .and. (CS%remap_uv_using_old_alg .and. (present(u) .or. present(v))) ) then - call MOM_error(FATAL, "remap_all_state_vars: dzInterface must be present if using old algorithm "// & - "and u/v are to be remapped") - endif - if (CS%answer_date >= 20190101) then h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff elseif (GV%Boussinesq) then @@ -876,7 +718,7 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 endif - if (show_call_tree) call callTree_enter("remap_all_state_vars(), MOM_ALE.F90") + if (show_call_tree) call callTree_enter("ALE_remap_tracers(), MOM_ALE.F90") nz = GV%ke @@ -890,7 +732,7 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & ! Remap all registered tracers, including temperature and salinity. if (ntr>0) then - if (show_call_tree) call callTree_waypoint("remapping tracers (remap_all_state_vars)") + if (show_call_tree) call callTree_waypoint("remapping tracers (ALE_remap_tracers)") !$OMP parallel do default(shared) private(h1,h2,tr_column,Tr,PCM,work_conc,work_cont,work_2d) do m=1,ntr ! For each tracer Tr => Reg%Tr(m) @@ -938,6 +780,7 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & if (Tr%id_remap_cont > 0) then call post_data(Tr%id_remap_cont, work_cont, CS%diag) endif + if (Tr%id_remap_cont_2d > 0) then do j = G%jsc,G%jec ; do i = G%isc,G%iec work_2d(i,j) = 0.0 @@ -952,9 +795,79 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & endif ! endif for ntr > 0 - if (show_call_tree) call callTree_waypoint("tracers remapped (remap_all_state_vars)") - if (CS%partial_cell_vel_remap .and. (present(u) .or. present(v)) ) then + if (CS%id_vert_remap_h > 0) call post_data(CS%id_vert_remap_h, h_old, CS%diag) + if ((CS%id_vert_remap_h_tendency > 0) .and. present(dt)) then + do k = 1, nz ; do j = G%jsc,G%jec ; do i = G%isc,G%iec + work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*Idt + enddo ; enddo ; enddo + call post_data(CS%id_vert_remap_h_tendency, work_cont, CS%diag) + endif + + if (show_call_tree) call callTree_leave("ALE_remap_tracers(), MOM_ALE.F90") + +end subroutine ALE_remap_tracers + +!> This routine remaps velocity components between the old and the new grids, +!! with thicknesses at velocity points taken to be arithmetic averages of tracer thicknesses. +!! This routine may be called during initialization of the model at time=0, to +!! remap initial conditions to the model grid. It is also called during a +!! time step to update the state. +subroutine ALE_remap_velocities(CS, G, GV, h_old, h_new, u, v, OBC, dzInterface, debug, dt) + type(ALE_CS), intent(in) :: CS !< ALE control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid + !! [H ~> m or kg m-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1] + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & + optional, intent(in) :: dzInterface !< Change in interface position + !! [H ~> m or kg m-2] + logical, optional, intent(in) :: debug !< If true, show the call tree + real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s] + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! The vertically summed thicknesses [H ~> m or kg m-2] + real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2] + real, dimension(GV%ke+1) :: dz ! The change in interface heights interpolated to + ! a velocity point [H ~> m or kg m-2] + logical :: PCM(GV%ke) ! If true, do PCM remapping from a cell. + real :: u_src(GV%ke) ! A column of u-velocities on the source grid [L T-1 ~> m s-1] + real :: u_tgt(GV%ke) ! A column of u-velocities on the target grid [L T-1 ~> m s-1] + real :: v_src(GV%ke) ! A column of v-velocities on the source grid [L T-1 ~> m s-1] + real :: v_tgt(GV%ke) ! A column of v-velocities on the target grid [L T-1 ~> m s-1] + real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] + real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2] + real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] + logical :: show_call_tree + integer :: i, j, k, nz + + show_call_tree = .false. + if (present(debug)) show_call_tree = debug + if (show_call_tree) call callTree_enter("ALE_remap_velocities()") + + ! If remap_uv_using_old_alg is .true. and u or v is requested, then we must have dzInterface. Otherwise, + ! u and v can be remapped without dzInterface + if (CS%remap_uv_using_old_alg .and. .not.present(dzInterface) ) call MOM_error(FATAL, & + "ALE_remap_velocities: dzInterface must be present if using old algorithm.") + + if (CS%answer_date >= 20190101) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + else + h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + endif + + nz = GV%ke + + if (CS%partial_cell_vel_remap) then h_tot(:,:) = 0.0 do k=1,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 h_tot(i,j) = h_tot(i,j) + h_old(i,j,k) @@ -962,13 +875,12 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & endif ! Remap u velocity component - if ( present(u) ) then + if ( .true. ) then !$OMP parallel do default(shared) private(h1,h2,dz,u_src,h_mask_vel,u_tgt) do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (G%mask2dCu(I,j)>0.) then ! Build the start and final grids do k=1,nz - u_src(k) = u(I,j,k) h1(k) = 0.5*(h_old(i,j,k) + h_old(i+1,j,k)) h2(k) = 0.5*(h_new(i,j,k) + h_new(i+1,j,k)) enddo @@ -994,6 +906,9 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & endif ; endif ! --- Remap u profiles from the source vertical grid onto the new target grid. + do k=1,nz + u_src(k) = u(I,j,k) + enddo call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt, & h_neglect, h_neglect_edge) @@ -1007,15 +922,14 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & endif ; enddo ; enddo endif - if (show_call_tree) call callTree_waypoint("u remapped (remap_all_state_vars)") + if (show_call_tree) call callTree_waypoint("u remapped (ALE_remap_velocities)") ! Remap v velocity component - if ( present(v) ) then + if ( .true. ) then !$OMP parallel do default(shared) private(h1,h2,v_src,dz,h_mask_vel,v_tgt) do J=G%JscB,G%JecB ; do i=G%isc,G%iec ; if (G%mask2dCv(i,J)>0.) then ! Build the start and final grids do k=1,nz - v_src(k) = v(i,J,k) h1(k) = 0.5*(h_old(i,j,k) + h_old(i,j+1,k)) h2(k) = 0.5*(h_new(i,j,k) + h_new(i,j+1,k)) enddo @@ -1039,6 +953,9 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & endif ; endif ! --- Remap v profiles from the source vertical grid onto the new target grid. + do k=1,nz + v_src(k) = v(i,J,k) + enddo call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt, & h_neglect, h_neglect_edge) @@ -1052,20 +969,13 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & endif ; enddo ; enddo endif - if (CS%id_vert_remap_h > 0) call post_data(CS%id_vert_remap_h, h_old, CS%diag) - if ((CS%id_vert_remap_h_tendency > 0) .and. present(dt)) then - do k = 1, nz ; do j = G%jsc,G%jec ; do i = G%isc,G%iec - work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*Idt - enddo ; enddo ; enddo - call post_data(CS%id_vert_remap_h_tendency, work_cont, CS%diag) - endif - if (show_call_tree) call callTree_waypoint("v remapped (remap_all_state_vars)") - if (show_call_tree) call callTree_leave("remap_all_state_vars()") + if (show_call_tree) call callTree_waypoint("v remapped (ALE_remap_velocities)") + if (show_call_tree) call callTree_leave("ALE_remap_velocities()") -end subroutine remap_all_state_vars +end subroutine ALE_remap_velocities -!> Mask out thicknesses to 0 when their runing sum exceeds a specified value. +!> Mask out thicknesses to 0 when their running sum exceeds a specified value. subroutine apply_partial_cell_mask(h1, h_mask) real, dimension(:), intent(inout) :: h1 !< A column of thicknesses to be masked out after their !! running vertical sum exceeds h_mask [H ~> m or kg m-2] @@ -1125,10 +1035,11 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c integer, intent(in) :: nk_src !< Number of levels on source grid real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: h_src !< Level thickness of source grid !! [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: s_src !< Scalar on source grid + real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: s_src !< Scalar on source grid, in arbitrary units [A] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid !! [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst !< Scalar on destination grid + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst !< Scalar on destination grid, in the same + !! arbitrary units as s_src [A] logical, optional, intent(in) :: all_cells !< If false, only reconstruct for !! non-vanished cells. Use all vanished !! layers otherwise (default). @@ -1142,8 +1053,8 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c !! for remapping ! Local variables integer :: i, j, k, n_points - real :: dx(GV%ke+1) - real :: h_neglect, h_neglect_edge + real :: dx(GV%ke+1) ! Change in interface position [H ~> m or kg m-2] + real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2] logical :: ignore_vanished_layers, use_remapping_core_w, use_2018_remap ignore_vanished_layers = .false. @@ -1222,18 +1133,18 @@ subroutine ALE_PLM_edge_values( CS, G, GV, h, Q, bdry_extrap, Q_t, Q_b ) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< layer thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: Q !< 3d scalar array + intent(in) :: Q !< 3d scalar array, in arbitrary units [A] logical, intent(in) :: bdry_extrap !< If true, use high-order boundary !! extrapolation within boundary cells real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: Q_t !< Scalar at the top edge of each layer + intent(inout) :: Q_t !< Scalar at the top edge of each layer [A] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: Q_b !< Scalar at the bottom edge of each layer + intent(inout) :: Q_b !< Scalar at the bottom edge of each layer [A] ! Local variables integer :: i, j, k - real :: slp(GV%ke) - real :: mslp - real :: h_neglect + real :: slp(GV%ke) ! Tracer slope times the cell width [A] + real :: mslp ! Monotonized tracer slope times the cell width [A] + real :: h_neglect ! Tiny thicknesses used in remapping [H ~> m or kg m-2] if (CS%answer_date >= 20190101) then h_neglect = GV%H_subroundoff @@ -1281,13 +1192,13 @@ subroutine TS_PPM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(ALE_CS), intent(inout) :: CS !< module control structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: S_t !< Salinity at the top edge of each layer + intent(inout) :: S_t !< Salinity at the top edge of each layer [S ~> ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: S_b !< Salinity at the bottom edge of each layer + intent(inout) :: S_b !< Salinity at the bottom edge of each layer [S ~> ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: T_t !< Temperature at the top edge of each layer + intent(inout) :: T_t !< Temperature at the top edge of each layer [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: T_b !< Temperature at the bottom edge of each layer + intent(inout) :: T_b !< Temperature at the bottom edge of each layer [C ~> degC] type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< layer thicknesses [H ~> m or kg m-2] @@ -1386,7 +1297,9 @@ end subroutine ALE_initRegridding function ALE_getCoordinate( CS ) type(ALE_CS), pointer :: CS !< module control structure - real, dimension(CS%nk+1) :: ALE_getCoordinate + real, dimension(CS%nk+1) :: ALE_getCoordinate !< The coordinate positions, in the appropriate units + !! of the target coordinate, e.g. [Z ~> m] for z*, + !! non-dimensional for sigma, etc. ALE_getCoordinate(:) = getCoordinateInterfaces( CS%regridCS, undo_scaling=.true. ) end function ALE_getCoordinate @@ -1416,7 +1329,7 @@ subroutine ALE_update_regrid_weights( dt, CS ) real, intent(in) :: dt !< Time-step used between ALE calls [T ~> s] type(ALE_CS), pointer :: CS !< ALE control structure ! Local variables - real :: w ! An implicit weighting estimate. + real :: w ! An implicit weighting estimate [nondim] if (associated(CS)) then w = 0.0 @@ -1443,7 +1356,7 @@ subroutine ALE_updateVerticalGridType(CS, GV) GV%zAxisUnits = getCoordinateUnits( CS%regridCS ) GV%zAxisLongName = getCoordinateShortName( CS%regridCS ) GV%direction = -1 ! Because of ferret in z* mode. Need method to set - ! as function of coordinae mode. + ! as function of coordinate mode. end subroutine ALE_updateVerticalGridType diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index e5ce4019ba..de287af98a 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -137,7 +137,7 @@ module MOM_regridding ! The following routines are visible to the outside world public initialize_regridding, end_regridding, regridding_main public regridding_preadjust_reqs, convective_adjustment -public inflate_vanished_layers_old, check_remapping_grid, check_grid_column +public inflate_vanished_layers_old, check_grid_column public set_regrid_params, get_regrid_size, write_regrid_file public uniformResolution, setCoordinateResolution public set_target_densities_from_GV, set_target_densities @@ -794,7 +794,7 @@ end subroutine end_regridding !------------------------------------------------------------------------------ !> Dispatching regridding routine for orchestrating regridding & remapping -subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, conv_adjust, & +subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, & frac_shelf_h, PCM_cell) !------------------------------------------------------------------------------ ! This routine takes care of (1) building a new grid and (2) remapping between @@ -823,24 +823,14 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, conv_ type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamical variables (T, S, ...) real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface - logical, intent(in ) :: conv_adjust !< If true, regridding_main should do - !! convective adjustment, but because it no - !! longer does convective adjustment this must - !! be false. This argument has been retained to - !! trap inconsistent code, but will eventually - !! be eliminated. real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out ) :: PCM_cell !< Use PCM remapping in cells where true ! Local variables real :: trickGnuCompiler + integer :: i, j - if (conv_adjust) call MOM_error(FATAL, & - "regridding_main: convective adjustment no longer is done inside of regridding_main. "//& - "The code needs to be modified to call regridding_main() with conv_adjust=.false, "//& - "and a call to convective_adjustment added before calling regridding_main() "//& - "if regridding_preadjust_reqs() indicates that this is necessary.") if (present(PCM_cell)) PCM_cell(:,:,:) = .false. select case ( CS%regridding_scheme ) @@ -879,8 +869,18 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, conv_ end select ! type of grid #ifdef __DO_SAFETY_CHECKS__ - if (CS%nk == GV%ke) call check_remapping_grid(G, GV, h, dzInterface,'in regridding_main') + if (CS%nk == GV%ke) then + do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 ; if (G%mask2dT(i,j)>0.) then + call check_grid_column( GV%ke, h(i,j,:), dzInterface(i,j,:), 'in regridding_main') + endif ; enddo ; enddo + endif #endif + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.) then + if (minval(h(i,j,:)) < 0.0) then + write(0,*) 'regridding_main check_grid: i,j=', i, j, 'h_new(i,j,:)=', h_new(i,j,:) + call MOM_error(FATAL, "regridding_main: negative thickness encountered.") + endif + endif ; enddo ; enddo end subroutine regridding_main @@ -952,23 +952,6 @@ subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) end subroutine calc_h_new_by_dz -!> Check that the total thickness of two grids match -subroutine check_remapping_grid( G, GV, h, dzInterface, msg ) - type(ocean_grid_type), intent(in) :: G !< Grid structure - type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzInterface !< Change in interface positions - !! [H ~> m or kg m-2] - character(len=*), intent(in) :: msg !< Message to append to errors - ! Local variables - integer :: i, j - - !$OMP parallel do default(shared) - do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 - if (G%mask2dT(i,j)>0.) call check_grid_column( GV%ke, h(i,j,:), dzInterface(i,j,:), msg ) - enddo ; enddo - -end subroutine check_remapping_grid !> Check that the total thickness of new and old grids are consistent subroutine check_grid_column( nk, h, dzInterface, msg ) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 267a162b00..c61f130ef7 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -50,9 +50,11 @@ module MOM use MOM_unit_tests, only : unit_tests ! MOM core modules -use MOM_ALE, only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity +use MOM_ALE, only : ALE_init, ALE_end, ALE_regrid, ALE_CS, adjustGridForIntegrity use MOM_ALE, only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile -use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags +use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, pre_ALE_adjustments +use MOM_ALE, only : ALE_remap_tracers, ALE_remap_velocities +use MOM_ALE, only : ALE_update_regrid_weights, pre_ALE_diagnostics, ALE_register_diags use MOM_ALE_sponge, only : rotate_ALE_sponge, update_ALE_sponge_field use MOM_barotropic, only : Barotropic_CS use MOM_boundary_update, only : call_OBC_register, OBC_register_end, update_OBC_CS @@ -160,7 +162,6 @@ module MOM use MOM_offline_main, only : offline_redistribute_residual, offline_diabatic_ale use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean use MOM_offline_main, only : offline_advection_layer, offline_transport_end -use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf use MOM_particles_mod, only : particles, particles_init, particles_run, particles_save_restart, particles_end @@ -1387,7 +1388,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) end subroutine step_MOM_tracer_dyn !> MOM_step_thermo orchestrates the thermodynamic time stepping and vertical -!! remapping, via calls to diabatic (or adiabatic) and ALE_main. +!! remapping, via calls to diabatic (or adiabatic) and ALE_regrid. subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & Time_end_thermo, update_BBL, Waves) type(MOM_control_struct), intent(inout) :: CS !< Master MOM control structure @@ -1409,13 +1410,17 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & optional, pointer :: Waves !< Container for wave related parameters !! the fields in Waves are intent in here. + real :: h_new(SZI_(G),SZJ_(G),SZK_(GV)) ! Layer thicknesses after regridding [H ~> m or kg m-2] + real :: dzRegrid(SZI_(G),SZJ_(G),SZK_(GV)+1) ! The change in grid interface positions due to regridding, + ! in the same units as thicknesses [H ~> m or kg m-2] + logical :: PCM_cell(SZI_(G),SZJ_(G),SZK_(GV)) ! If true, PCM remapping should be used in a cell. logical :: use_ice_shelf ! Needed for selecting the right ALE interface. logical :: showCallTree type(group_pass_type) :: pass_T_S, pass_T_S_h, pass_uv_T_S_h integer :: dynamics_stencil ! The computational stencil for the calculations ! in the dynamic core. integer :: halo_sz ! The size of a halo where data must be valid. - integer :: is, ie, js, je, nz + integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke showCallTree = callTree_showQuery() @@ -1486,7 +1491,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! Regridding/remapping is done here, at end of thermodynamics time step ! (that may comprise several dynamical time steps) - ! The routine 'ALE_main' can be found in 'MOM_ALE.F90'. + ! The routine 'ALE_regrid' can be found in 'MOM_ALE.F90'. if ( CS%use_ALE_algorithm ) then call enable_averages(dtdia, Time_end_thermo, CS%diag) ! call pass_vector(u, v, G%Domain) @@ -1508,14 +1513,32 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call check_redundant("Pre-ALE ", u, v, G) endif call cpu_clock_begin(id_clock_ALE) + + call pre_ALE_diagnostics(G, GV, US, h, u, v, tv, CS%ALE_CSp) + call ALE_update_regrid_weights(dtdia, CS%ALE_CSp) + ! Do any necessary adjustments ot the state prior to remapping. + call pre_ALE_adjustments(G, GV, US, h, tv, CS%tracer_Reg, CS%ALE_CSp, u, v) + ! Adjust the target grids for diagnostics, in case there have been thickness adjustments. + call diag_update_remap_grids(CS%diag) + if (use_ice_shelf) then - call ALE_main(G, GV, US, h, u, v, tv, CS%tracer_Reg, CS%ALE_CSp, CS%OBC, & - dtdia, CS%frac_shelf_h) + call ALE_regrid(G, GV, US, h, h_new, dzRegrid, tv, CS%ALE_CSp, CS%frac_shelf_h, PCM_cell) else - call ALE_main(G, GV, US, h, u, v, tv, CS%tracer_Reg, CS%ALE_CSp, CS%OBC, dtdia) + call ALE_regrid(G, GV, US, h, h_new, dzRegrid, tv, CS%ALE_CSp, PCM_cell=PCM_cell) endif - if (showCallTree) call callTree_waypoint("finished ALE_main (step_MOM_thermo)") + if (showCallTree) call callTree_waypoint("new grid generated") + ! Remap all variables from the old grid h onto the new grid h_new + call ALE_remap_tracers(CS%ALE_CSp, G, GV, h, h_new, CS%tracer_Reg, showCallTree, dtdia, PCM_cell) + call ALE_remap_velocities(CS%ALE_CSp, G, GV, h, h_new, u, v, CS%OBC, dzRegrid, showCallTree, dtdia) + + ! Replace the old grid with new one. All remapping must be done by this point in the code. + !$OMP parallel do default(shared) + do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + h(i,j,k) = h_new(i,j,k) + enddo ; enddo ; enddo + + if (showCallTree) call callTree_waypoint("finished ALE_regrid (step_MOM_thermo)") call cpu_clock_end(id_clock_ALE) endif ! endif for the block "if ( CS%use_ALE_algorithm )" @@ -1614,13 +1637,16 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS logical :: do_vertical !< If enough time has elapsed, do the diabatic tracer sources/sinks logical :: adv_converged !< True if all the horizontal fluxes have been used + real, allocatable, dimension(:,:,:) :: h_new ! Layer thicknesses after regridding [H ~> m or kg m-2] + real, allocatable, dimension(:,:,:) :: dzRegrid ! The change in grid interface positions due to regridding, + ! in the same units as thicknesses [H ~> m or kg m-2] real :: dt_offline ! The offline timestep for advection [T ~> s] real :: dt_offline_vertical ! The offline timestep for vertical fluxes and remapping [T ~> s] logical :: skip_diffusion type(time_type), pointer :: accumulated_time => NULL() type(time_type), pointer :: vertical_time => NULL() - integer :: is, ie, js, je, isd, ied, jsd, jed + integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz ! 3D pointers real, dimension(:,:,:), pointer :: & @@ -1635,7 +1661,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS ! Grid-related pointer assignments G => CS%G ; GV => CS%GV ; US => CS%US - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed call cpu_clock_begin(id_clock_offline_tracer) @@ -1646,19 +1672,11 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call enable_averaging(time_interval, Time_end, CS%diag) ! Check to see if this is the first iteration of the offline interval - if (accumulated_time == real_to_time(0.0)) then - first_iter = .true. - else ! This is probably unnecessary but is used to guard against unwanted behavior - first_iter = .false. - endif + first_iter = (accumulated_time == real_to_time(0.0)) ! Check to see if vertical tracer functions should be done - if (first_iter .or. (accumulated_time >= vertical_time)) then - do_vertical = .true. - vertical_time = accumulated_time + real_to_time(US%T_to_s*dt_offline_vertical) - else - do_vertical = .false. - endif + do_vertical = (first_iter .or. (accumulated_time >= vertical_time)) + if (do_vertical) vertical_time = accumulated_time + real_to_time(US%T_to_s*dt_offline_vertical) ! Increment the amount of time elapsed since last read and check if it's time to roll around accumulated_time = accumulated_time + real_to_time(time_interval) @@ -1737,7 +1755,28 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS ! Call ALE one last time to make sure that tracers are remapped onto the layer thicknesses ! stored from the forward run call cpu_clock_begin(id_clock_ALE) - call ALE_offline_tracer_final( G, GV, CS%h, CS%tv, h_end, CS%tracer_Reg, CS%ALE_CSp, CS%OBC) + + ! Do any necessary adjustments ot the state prior to remapping. + call pre_ALE_adjustments(G, GV, US, h_end, CS%tv, CS%tracer_Reg, CS%ALE_CSp) + + allocate(h_new(isd:ied, jsd:jed, nz), source=0.0) + allocate(dzRegrid(isd:ied, jsd:jed, nz+1), source=0.0) + + ! Generate the new grid based on the tracer grid at the end of the interval. + call ALE_regrid(G, GV, US, h_end, h_new, dzRegrid, CS%tv, CS%ALE_CSp) + + ! Remap the tracers from the previous tracer grid onto the new grid. The thicknesses that + ! are used are intended to ensure that in the case where transports don't quite conserve, + ! the offline layer thicknesses do not drift too far away from the online model. + call ALE_remap_tracers(CS%ALE_CSp, G, GV, CS%h, h_new, CS%tracer_Reg, debug=CS%debug) + + ! Update the tracer grid. + do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + CS%h(i,j,k) = h_new(i,j,k) + enddo ; enddo ; enddo + + deallocate(h_new, dzRegrid) + call cpu_clock_end(id_clock_ALE) call pass_var(CS%h, G%Domain) endif @@ -1847,12 +1886,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! This include declares and sets the variable "version". # include "version_variable.h" - integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB real :: dtbt ! If negative, this specifies the barotropic timestep as a fraction ! of the maximum stable value [nondim]. real, allocatable, dimension(:,:) :: eta ! free surface height or column mass [H ~> m or kg m-2] + real, allocatable, dimension(:,:,:) :: h_new ! Layer thicknesses after regridding [H ~> m or kg m-2] + real, allocatable, dimension(:,:,:) :: dzRegrid ! The change in grid interface positions due to regridding, + ! in the same units as thicknesses [H ~> m or kg m-2] + logical, allocatable, dimension(:,:,:) :: PCM_cell ! If true, PCM remapping should be used in a cell. type(group_pass_type) :: tmp_pass_uv_T_S_h, pass_uv_T_S_h real :: default_val ! default value for a parameter @@ -2775,14 +2818,30 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif call callTree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)") call adjustGridForIntegrity(CS%ALE_CSp, G, GV, CS%h ) - call callTree_waypoint("Calling ALE_main() to remap initial conditions (initialize_MOM)") + call pre_ALE_adjustments(G, GV, US, CS%h, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%u, CS%v) + + call callTree_waypoint("Calling ALE_regrid() to remap initial conditions (initialize_MOM)") + allocate(h_new(isd:ied, jsd:jed, nz), source=0.0) + allocate(dzRegrid(isd:ied, jsd:jed, nz+1), source=0.0) + allocate(PCM_cell(isd:ied, jsd:jed, nz), source=.false.) if (use_ice_shelf) then - call ALE_main(G, GV, US, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, & - CS%OBC, frac_shelf_h=CS%frac_shelf_h) + call ALE_regrid(G, GV, US, CS%h, h_new, dzRegrid, CS%tv, CS%ALE_CSp, CS%frac_shelf_h, PCM_cell) else - call ALE_main( G, GV, US, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%OBC) + call ALE_regrid(G, GV, US, CS%h, h_new, dzRegrid, CS%tv, CS%ALE_CSp, PCM_cell=PCM_cell) endif + if (callTree_showQuery()) call callTree_waypoint("new grid generated") + ! Remap all variables from the old grid h onto the new grid h_new + call ALE_remap_tracers(CS%ALE_CSp, G, GV, CS%h, h_new, CS%tracer_Reg, CS%debug, PCM_cell=PCM_cell) + call ALE_remap_velocities(CS%ALE_CSp, G, GV, CS%h, h_new, CS%u, CS%v, CS%OBC, dzRegrid, debug=CS%debug) + + ! Replace the old grid with new one. All remapping must be done at this point. + !$OMP parallel do default(shared) + do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + CS%h(i,j,k) = h_new(i,j,k) + enddo ; enddo ; enddo + deallocate(h_new, dzRegrid, PCM_cell) + call cpu_clock_begin(id_clock_pass_init) call create_group_pass(tmp_pass_uv_T_S_h, CS%u, CS%v, G%Domain) if (use_temperature) then @@ -2796,6 +2855,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (CS%debug) then call uvchksum("Post ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1) call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_m) + if (use_temperature) then + call hchksum(CS%tv%T, "Post ALE adjust init cond T", G%HI, haloshift=1, scale=US%C_to_degC) + call hchksum(CS%tv%S, "Post ALE adjust init cond S", G%HI, haloshift=1, scale=US%S_to_ppt) + endif endif endif if ( CS%use_ALE_algorithm ) call ALE_updateVerticalGridType( CS%ALE_CSp, GV ) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0eac15f6d7..31dbb41dcc 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2779,7 +2779,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just call regridding_preadjust_reqs(regridCS, do_conv_adj, ignore) if (do_conv_adj) call convective_adjustment(G, GV_loc, h1, tv_loc) - call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, conv_adjust=.false., & + call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, & frac_shelf_h=frac_shelf_h ) deallocate( dz_interface ) diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index c1582dca4a..bf06fc294e 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -4,7 +4,9 @@ module MOM_offline_main ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_ALE, only : ALE_CS, ALE_main_offline, ALE_offline_inputs +use MOM_ALE, only : ALE_CS, ALE_regrid, ALE_offline_inputs +use MOM_ALE, only : pre_ALE_adjustments, ALE_update_regrid_weights +use MOM_ALE, only : ALE_remap_tracers use MOM_checksums, only : hchksum, uvchksum use MOM_coms, only : reproducing_sum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end @@ -118,7 +120,7 @@ module MOM_offline_main real :: minimum_forcing_depth !< The smallest depth over which fluxes can be applied [H ~> m or kg m-2]. !! This is copied from diabatic_CS controlling how tracers follow freshwater fluxes - real :: Kd_max !< Runtime parameter specifying the maximum value of vertical diffusivity + real :: Kd_max !< Runtime parameter specifying the maximum value of vertical diffusivity [Z2 T-1 ~> m2 s-1] real :: min_residual !< The minimum amount of total mass flux before exiting the main advection !! routine [H L2 ~> m3 or kg] !>@{ Diagnostic manager IDs for some fields that may be of interest when doing offline transport @@ -169,8 +171,6 @@ module MOM_offline_main real, allocatable, dimension(:,:,:) :: Kd !< Vertical diffusivity [Z2 T-1 ~> m2 s-1] real, allocatable, dimension(:,:,:) :: h_end !< Thicknesses at the end of offline timestep [H ~> m or kg m-2] - real, allocatable, dimension(:,:) :: netMassIn !< Freshwater fluxes into the ocean - real, allocatable, dimension(:,:) :: netMassOut !< Freshwater fluxes out of the ocean real, allocatable, dimension(:,:) :: mld !< Mixed layer depths at thickness points [Z ~> m] ! Allocatable arrays to read in entire fields during initialization @@ -229,7 +229,10 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C ! Variables used to keep track of layer thicknesses at various points in the code real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & h_new, & ! Updated layer thicknesses [H ~> m or kg m-2] + h_post_remap, & ! Layer thicknesses after remapping [H ~> m or kg m-2] h_vol ! Layer volumes [H L2 ~> m3 or kg] + real :: dzRegrid(SZI_(G),SZJ_(G),SZK_(GV)+1) ! The change in grid interface positions due to regridding, + ! in the same units as thicknesses [H ~> m or kg m-2] integer :: niter, iter real :: Inum_iter ! The inverse of the number of iterations [nondim] character(len=256) :: mesg ! The text of an error message @@ -347,7 +350,23 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C call MOM_tracer_chkinv(debug_msg, G, GV, h_new, CS%tracer_reg) endif call cpu_clock_begin(id_clock_ALE) - call ALE_main_offline(G, GV, h_new, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%OBC, CS%dt_offline) + + call ALE_update_regrid_weights(CS%dt_offline, CS%ALE_CSp) + call pre_ALE_adjustments(G, GV, US, h_new, CS%tv, CS%tracer_Reg, CS%ALE_CSp) + ! Uncomment this to adjust the target grids for diagnostics, if there have been thickness + ! adjustments, but the offline tracer code does not yet have the other corresponding calls + ! that would be needed to support remapping its output. + ! call diag_update_remap_grids(CS%diag, alt_h=h_new) + + call ALE_regrid(G, GV, US, h_new, h_post_remap, dzRegrid, CS%tv, CS%ALE_CSp) + + ! Remap all variables from the old grid h_new onto the new grid h_post_remap + call ALE_remap_tracers(CS%ALE_CSp, G, GV, h_new, h_post_remap, CS%tracer_Reg, & + CS%debug, dt=CS%dt_offline) + + do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + h_new(i,j,k) = h_post_remap(i,j,k) + enddo ; enddo ; enddo call cpu_clock_end(id_clock_ALE) if (CS%debug) then @@ -746,6 +765,7 @@ subroutine offline_fw_fluxes_into_ocean(G, GV, CS, fluxes, h, in_flux_optional) real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: in_flux_optional !< The total time-integrated amount !! of tracer that leaves with freshwater + !! [CU H ~> Conc m or Conc kg m-2] integer :: i, j, m real, dimension(SZI_(G),SZJ_(G)) :: negative_fw !< store all negative fluxes [H ~> m or kg m-2] @@ -796,6 +816,7 @@ subroutine offline_fw_fluxes_out_ocean(G, GV, CS, fluxes, h, out_flux_optional) real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: out_flux_optional !< The total time-integrated amount !! of tracer that leaves with freshwater + !! [CU H ~> Conc m or Conc kg m-2] integer :: m logical :: update_h !< Flag for whether h should be updated @@ -1444,8 +1465,6 @@ subroutine offline_transport_init(param_file, CS, diabatic_CSp, G, GV, US) allocate(CS%eatr(isd:ied,jsd:jed,nz), source=0.0) allocate(CS%ebtr(isd:ied,jsd:jed,nz), source=0.0) allocate(CS%h_end(isd:ied,jsd:jed,nz), source=0.0) - allocate(CS%netMassOut(G%isd:G%ied,G%jsd:G%jed), source=0.0) - allocate(CS%netMassIn(G%isd:G%ied,G%jsd:G%jed), source=0.0) allocate(CS%Kd(isd:ied,jsd:jed,nz+1), source=0.0) if (CS%read_mld) allocate(CS%mld(G%isd:G%ied,G%jsd:G%jed), source=0.0) @@ -1518,8 +1537,6 @@ subroutine offline_transport_end(CS) deallocate(CS%eatr) deallocate(CS%ebtr) deallocate(CS%h_end) - deallocate(CS%netMassOut) - deallocate(CS%netMassIn) deallocate(CS%Kd) if (CS%read_mld) deallocate(CS%mld) if (CS%read_all_ts_uvh) then