Merge branch 'ashao-offline_tracer_updates' into dev/gfdl

config_src/coupled_driver/ocean_model_MOM.F90

@@ -36,7 +36,7 @@ module ocean_model_mod
 
 use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end
 use MOM, only : calculate_surface_state, finish_MOM_initialization
-use MOM, only : step_tracers
+use MOM, only : step_offline
 use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf
 use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging
 use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end
@@ -470,7 +470,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
   Master_time = OS%Time ; Time1 = OS%Time
 
   if(OS%MOM_Csp%offline_tracer_mode) then
-    call step_tracers(OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp)
+    call step_offline(OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp)
   else
     call step_MOM(OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp)
   endif

config_src/solo_driver/MOM_driver.F90

@@ -47,7 +47,7 @@ program MOM_main
   use MOM_diag_mediator,   only : diag_mediator_close_registration, diag_mediator_end
   use MOM,                 only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end
   use MOM,                 only : calculate_surface_state, finish_MOM_initialization
-  use MOM,                 only : step_tracers
+  use MOM,                 only : step_offline
   use MOM_domains,         only : MOM_infra_init, MOM_infra_end
   use MOM_error_handler,   only : MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe
   use MOM_error_handler,   only : callTree_enter, callTree_leave, callTree_waypoint
@@ -464,7 +464,7 @@ program MOM_main
     ! This call steps the model over a time time_step.
     Time1 = Master_Time ; Time = Master_Time
     if (offline_tracer_mode) then
-      call step_tracers(fluxes, state, Time1, time_step, MOM_CSp)
+      call step_offline(fluxes, state, Time1, time_step, MOM_CSp)
     else
       call step_MOM(fluxes, state, Time1, time_step, MOM_CSp)
     endif

src/ALE/MOM_ALE.F90

@@ -8,9 +8,12 @@ module MOM_ALE
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use MOM_debugging,        only : check_column_integrals
 use MOM_diag_mediator,    only : register_diag_field, post_data, diag_ctrl, time_type
+use MOM_diag_vkernels,    only : interpolate_column, reintegrate_column
 use MOM_domains,          only : create_group_pass, do_group_pass, group_pass_type
 use MOM_EOS,              only : calculate_density
+use MOM_domains,          only : create_group_pass, do_group_pass, group_pass_type
 use MOM_error_handler,    only : MOM_error, FATAL, WARNING
 use MOM_error_handler,    only : callTree_showQuery
 use MOM_error_handler,    only : callTree_enter, callTree_leave, callTree_waypoint
@@ -36,7 +39,7 @@ module MOM_ALE
 use MOM_remapping,        only : remappingSchemesDoc, remappingDefaultScheme
 use MOM_remapping,        only : remapping_CS, dzFromH1H2
 use MOM_string_functions, only : uppercase, extractWord, extract_integer
-use MOM_tracer_registry,  only : tracer_registry_type
+use MOM_tracer_registry,  only : tracer_registry_type, MOM_tracer_chkinv
 use MOM_variables,        only : ocean_grid_type, thermo_var_ptrs
 use MOM_verticalGrid,     only : verticalGrid_type
 
@@ -105,6 +108,7 @@ module MOM_ALE
 public ALE_end
 public ALE_main
 public ALE_main_offline
+public ALE_offline_inputs
 public ALE_offline_tracer_final
 public ALE_build_grid
 public ALE_regrid_accelerated
@@ -488,51 +492,130 @@ subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, dt)
 
 end subroutine ALE_main_offline
 
+!> 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
+!! routine builds a grid on the runtime specified vertical coordinate
+subroutine ALE_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug)
+  type(ALE_CS),                                 pointer       :: CS    !< Regridding parameters and options
+  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     !< Layer thicknesses
+  type(thermo_var_ptrs),                        intent(inout) :: tv    !< Thermodynamic variable structure
+  type(tracer_registry_type),                   pointer       :: Reg   !< Tracer registry structure
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)),   intent(inout) :: uhtr  !< Zonal mass fluxes
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),   intent(inout) :: vhtr  !< Meridional mass fluxes
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),  intent(inout) :: Kd    !< Input diffusivites
+  logical,                                      intent(in   ) :: debug !< If true, then turn checksums
+  ! 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
+  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions
+  real, dimension(SZK_(GV)) :: h_src
+  real, dimension(SZK_(GV)) :: h_dest, uh_dest
+  real, dimension(SZK_(GV)) :: temp_vec
+
+  nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
+  dzRegrid(:,:,:) = 0.0
+  h_new(:,:,:) = 0.0
+
+  if (debug) call MOM_tracer_chkinv("Before ALE_offline_inputs", G, h, Reg%Tr, Reg%ntr)
+
+  ! 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. )
+  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%remapCS, 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
+  do j=jsc,jec ; do i=G%iscB,G%iecB
+    if (G%mask2dCu(i,j)>0.) then
+      h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
+      h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i+1,j,:))
+      call reintegrate_column(nk, h_src, uhtr(I,j,:), nk, h_dest, 0., temp_vec)
+      uhtr(I,j,:) = temp_vec
+    endif
+  enddo ; enddo
+  do j=G%jscB,G%jecB ; do i=isc,iec
+    if (G%mask2dCv(i,j)>0.) then
+      h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
+      h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i,j+1,:))
+      call reintegrate_column(nk, h_src, vhtr(I,j,:), nk, h_dest, 0., temp_vec)
+      vhtr(I,j,:) = temp_vec
+    endif
+  enddo ; enddo
+
+  do j = jsc,jec ; do i=isc,iec
+    if (G%mask2dT(i,j)>0.) then
+      if (check_column_integrals(nk, h_src, nk, h_dest)) then
+        call MOM_error(FATAL, "ALE_offline_inputs: Kd interpolation columns do not match")
+      endif
+      call interpolate_column(nk, h(i,j,:), Kd(i,j,:), nk, h_new(i,j,:), 0., Kd(i,j,:))
+    endif
+  enddo ; enddo;
+
+  call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%T, h_new, tv%T)
+  call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%S, h_new, tv%S)
+
+  if (debug) call MOM_tracer_chkinv("After ALE_offline_inputs", G, h_new, Reg%Tr, Reg%ntr)
+
+  ! 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)
+  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, h_target, Reg, CS)
+subroutine ALE_offline_tracer_final( G, GV, h, tv, h_target, Reg, CS)
   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 last time step (m or Pa)
-  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(in)    :: h_target   !< Current 3D grid obtained after last time step (m or Pa)
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  intent(inout) :: h   !< Current 3D grid obtained after
+                                                                   !! last time step (m or Pa)
+  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 (m or Pa)
   type(tracer_registry_type),                 pointer       :: Reg !< Tracer registry structure
   type(ALE_CS),                               pointer       :: CS  !< Regridding parameters and options
   ! 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) :: dzRegrid !< The change in grid interface positions
+  real, dimension(SZI_(G), SZJ_(G), SZK_(GV))   :: h_new    !< Regridded target thicknesses
   integer :: nk, i, j, k, isc, iec, jsc, jec
 
   nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
 
   if (CS%show_call_tree) call callTree_enter("ALE_offline_tracer_final(), MOM_ALE.F90")
-
-  ! It does not seem that remap_all_state_vars uses dzRegrid for tracers, only for u, v
-  dzRegrid(:,:,:) = 0.0
-
-  call check_grid( G, GV, h, 0. )
+  ! Need to make sure that h_target is consistent with the current offline ALE confiuration
+  call regridding_main( CS%remapCS, CS%regridCS, G, GV, h_target, tv, h_new, dzRegrid )
   call check_grid( G, GV, h_target, 0. )
 
-  if (CS%show_call_tree) call callTree_waypoint("Source and target grids checked (ALE_offline_tracer)")
+
+  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%remapCS, CS, G, GV, h, h_target, Reg, &
-                             debug=CS%show_call_tree )
+  call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, Reg, debug=CS%show_call_tree )
 
-  if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_offline_tracer)")
+  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(none) shared(isc,iec,jsc,jec,nk,h,h_target,CS)
   do k = 1,nk
     do j = jsc-1,jec+1 ; do i = isc-1,iec+1
-      h(i,j,k) = h_target(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_tracer()")
-
+  if (CS%show_call_tree) call callTree_leave("ALE_offline_tracer_final()")
 end subroutine ALE_offline_tracer_final
 
 !> Check grid for negative thicknesses

src/ALE/MOM_regridding.F90

@@ -757,7 +757,7 @@ end subroutine end_regridding
 !------------------------------------------------------------------------------
 ! Dispatching regridding routine: regridding & remapping
 !------------------------------------------------------------------------------
-subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h)
+subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
 !------------------------------------------------------------------------------
 ! This routine takes care of (1) building a new grid and (2) remapping between
 ! the old grid and the new grid. The creation of the new grid can be based
@@ -785,9 +785,14 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_
   real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h_new  !< New 3D grid consistent with target coordinate
   real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in position of each interface
   real, dimension(:,:),                   optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
+  logical,                          optional, intent(in   ) :: conv_adjust ! If true, do convective adjustment
   ! Local variables
   real :: trickGnuCompiler
   logical :: use_ice_shelf
+  logical :: do_convective_adjustment
+
+  do_convective_adjustment = .true.
+  if (present(conv_adjust)) do_convective_adjustment = conv_adjust
 
   use_ice_shelf = .false.
   if (present(frac_shelf_h)) then
@@ -813,7 +818,7 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_
       call calc_h_new_by_dz(G, GV, h, dzInterface, h_new)
 
     case ( REGRIDDING_RHO )
-      call convective_adjustment(G, GV, h, tv)
+      if (do_convective_adjustment) call convective_adjustment(G, GV, h, tv)
       call build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
       call calc_h_new_by_dz(G, GV, h, dzInterface, h_new)
 
@@ -1275,6 +1280,9 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
   integer :: i, j, k
   real    :: nominalDepth, totalThickness
   real, dimension(SZK_(GV)+1) :: zOld, zNew
+#ifdef __DO_SAFETY_CHECKS__
+  real    :: dh
+#endif
 
   nz = GV%ke