Changes made so that regridding onto pure isopycnal kind of works. Ne…

…ed to

figure out why the redistribution of residual does not fully conserve.

src/ALE/MOM_ALE.F90

@@ -8,6 +8,7 @@ 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
@@ -510,6 +511,7 @@ subroutine ALE_offline_inputs(CS, G, GV, h_input, h_regrid, tv, Reg, uhtr, vhtr,
 
   nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec
   dzRegrid(:,:,:) = 0.0
+  h_regrid(:,:,:) = 0.0
 
   ! 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.
@@ -523,18 +525,33 @@ subroutine ALE_offline_inputs(CS, G, GV, h_input, h_regrid, tv, Reg, uhtr, vhtr,
 
   ! Reintegrate mass transports
   do j=jsc,jec ; do i=G%iscB,G%iecB
-    h_src(:) = 0.5 * (h_input(i,j,:) + h_input(i+1,j,:))
-    h_dest(:) = 0.5 * (h_regrid(i,j,:) + h_regrid(i+1,j,:))
-    call reintegrate_column(nk, h_src, uhtr(I,j,:), nk, h_dest, 0., uhtr(I,j,:))
+    if (G%mask2dCu(i,j)>0.) then
+      h_src(:) = 0.5 * (h_input(i,j,:) + h_input(i+1,j,:))
+      h_dest(:) = 0.5 * (h_regrid(i,j,:) + h_regrid(i+1,j,:))
+      if (check_column_integrals(nk, h_src, nk, h_dest)) then
+        call MOM_error(FATAL, "ALE_offline_inputs: uh reintegration columns do not match")
+      endif
+      call reintegrate_column(nk, h_src, uhtr(I,j,:), nk, h_dest, 0., uhtr(I,j,:))
+    endif
   enddo ; enddo
   do j=G%jscB,G%jecB ; do i=isc,iec
-    h_src(:) = 0.5 * (h_input(i,j,:) + h_input(i,j+1,:))
-    h_dest(:) = 0.5 * (h_regrid(i,j,:) + h_regrid(i,j+1,:))
-    call reintegrate_column(nk, h_src, vhtr(I,j,:), nk, h_dest, 0., vhtr(I,j,:))
+    if (G%mask2dCv(i,j)>0.) then
+      h_src(:) = 0.5 * (h_input(i,j,:) + h_input(i,j+1,:))
+      h_dest(:) = 0.5 * (h_regrid(i,j,:) + h_regrid(i,j+1,:))
+      if (check_column_integrals(nk, h_src, nk, h_dest)) then
+        call MOM_error(FATAL, "ALE_offline_inputs: vh reintegration columns do not match")
+      endif
+      call reintegrate_column(nk, h_src, vhtr(I,j,:), nk, h_dest, 0., vhtr(I,j,:))
+    endif
   enddo ; enddo
 
   do j = jsc,jec ; do i=isc,iec
-    call interpolate_column(nk, h_input(i,j,:), Kd(i,j,:), nk, h_regrid(i,j,:), 0., Kd(i,j,:))
+    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_input(i,j,:), Kd(i,j,:), nk, h_regrid(i,j,:), 0., Kd(i,j,:))
+    endif
   enddo ; enddo;
 
   call ALE_remap_scalar(CS%remapCS, G, GV, nk, h_input, tv%T, h_regrid, tv%T)

src/core/MOM.F90

@@ -1564,7 +1564,7 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS)
   ! 3D pointers
   real, dimension(:,:,:), pointer   :: &
     uhtr, vhtr, &
-    eatr, ebtr, &    
+    eatr, ebtr, &
     h_end
 
   ! 2D Array for diagnostics
@@ -1707,14 +1707,13 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS)
       call tracer_hordiff(h_end, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, &
         CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv)
     endif
-    
+
     CS%h = h_end
 
     call pass_var(CS%tv%T,G%Domain)
     call pass_var(CS%tv%S,G%Domain)
     call pass_var(CS%h,G%Domain)
 
-
   endif
 
   call adjust_ssh_for_p_atm(CS, G, GV, CS%ave_ssh, fluxes%p_surf_SSH)

src/tracer/MOM_offline_aux.F90

@@ -651,6 +651,16 @@ subroutine update_offline_from_files(G, GV, nk_input, mean_file, sum_file, snap_
       timelevel=ridx_sum,position=CENTER)
   endif
 
+  do j=js,je ; do i=is,ie
+    if (G%mask2dT(i,j)>0.) then
+      do k=1,nz
+        if (h_start(i,j,k)<GV%Angstrom) exit
+      enddo
+      temp_mean(:,:,k:nz) = temp_mean(i,j,k)
+      salt_mean(:,:,k:nz) = salt_mean(i,j,k)
+    endif
+  enddo ; enddo 
+
   ! Check if reading vertical diffusivities or entrainment fluxes
   call read_data( mean_file, 'Kd_interface', Kd(:,:,1:nk_input+1), domain=G%Domain%mpp_domain, &
                   timelevel=ridx_sum,position=CENTER)

src/tracer/MOM_offline_main.F90

@@ -116,6 +116,14 @@ module MOM_offline_main
     id_h_redist = -1, &
     id_eta_diff_end = -1
 
+  !> Diagnostic IDs for the regridded/remapped input fields  
+  integer :: &
+    id_uhtr_regrid = -1, &
+    id_vhtr_regrid = -1, &
+    id_temp_regrid = -1, &
+    id_salt_regrid = -1, &
+    id_h_regrid = -1
+
   !> IDs for timings of various offline components
   integer ::  &
     id_clock_read_fields = -1,      &
@@ -272,8 +280,12 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, CS, &
     call hchksum(h_pre,"h_pre before transport",G%HI)
     call uvchksum("[uv]htr_sub before transport", uhtr_sub, vhtr_sub, G%HI)
   endif
-  ! Initialize the previous tot_residual to a negative value; tot_residual will always be positive
-  prev_tot_residual = -1.
+  tot_residual = remaining_transport_sum(CS, uhtr, vhtr)
+  if (CS%print_adv_offline) then
+    if (is_root_pe()) then
+      write(*,'(A,ES24.16)') "Main advection starting transport: ", tot_residual
+    endif
+  endif
 
   ! This loop does essentially a flux-limited, nonlinear advection scheme until all mass fluxes
   ! are used. ALE is done after the horizontal advection.
@@ -305,7 +317,7 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, CS, &
       h_new(i,j,k) = h_new(i,j,k)/G%areaT(i,j)
     enddo ; enddo ; enddo
 
-    if (MODULO(iter,1)==0) then
+    if (MODULO(iter,100)==0) then
       ! Do ALE remapping/regridding to allow for more advection to occur in the next iteration
       call pass_var(h_new,G%Domain)
       if (CS%debug) then
@@ -399,7 +411,7 @@ subroutine offline_redistribute_residual(CS, h_pre, uhtr, vhtr, converged, id_cl
   real, dimension(SZI_(CS%G),SZJB_(CS%G),SZK_(CS%G)) :: vhr  !< Meridional mass transport
 
   integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz, iter
-  real :: tot_residual, stock_values(MAX_FIELDS_)
+  real :: prev_tot_residual, tot_residual, stock_values(MAX_FIELDS_)
   integer :: nstocks
 
   ! Assign grid pointers
@@ -535,9 +547,13 @@ subroutine offline_redistribute_residual(CS, h_pre, uhtr, vhtr, converged, id_cl
         endif
       endif
       ! If the remaining residual is 0, then this return is done
-      if (tot_residual==0.0) then
+      if (tot_residual==0.0 ) then
         exit
       endif
+
+      if ( (tot_residual == prev_tot_residual) .or. (tot_residual<CS%min_residual) ) exit
+      prev_tot_residual = tot_residual
+
     enddo ! Redistribution iteration
   endif ! If one of the redistribution routines is requested
 
@@ -949,6 +965,7 @@ subroutine update_offline_fields(CS, h, fluxes, do_ale)
   logical,              intent(in   ) :: do_ale !< True if using ALE
   ! Local variables
   integer :: i, j, k, is, ie, js, je, nz
+  real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: temp_mean
   is = CS%G%isc ; ie = CS%G%iec ; js = CS%G%jsc ; je = CS%G%jec ; nz = CS%GV%ke
 
   call cpu_clock_begin(CS%id_clock_read_fields)
@@ -961,14 +978,26 @@ subroutine update_offline_fields(CS, h, fluxes, do_ale)
   call update_offline_from_arrays(CS%G, CS%GV, CS%nk_input, CS%ridx_sum, CS%mean_file, CS%sum_file, CS%snap_file,   &
                                   CS%uhtr, CS%vhtr, CS%h_start, CS%h_end, CS%uhtr_all, CS%vhtr_all, CS%hstart_all,  &
                                   CS%hend_all, CS%tv%T, CS%tv%s, CS%temp_all, CS%salt_all, CS%read_all_ts_uvh)
-  call uvchksum("[uv]h after update offline from files and arrays", CS%uhtr, CS%vhtr, CS%G%HI)
-
+  if (CS%debug) then
+    call uvchksum("[uv]h after update offline from files and arrays", CS%uhtr, CS%vhtr, CS%G%HI)  
+  endif
+
   ! These halo passes are necessary because u, v fields will need information 1 step into the halo
   call pass_var(CS%h_start, CS%G%Domain)
   call pass_var(CS%tv%T, CS%G%Domain)
   call pass_var(CS%tv%S, CS%G%Domain)
   call ALE_offline_inputs(CS%ALE_CSp, CS%G, CS%GV,  CS%h_start, h, CS%tv, CS%tracer_Reg, CS%uhtr, CS%vhtr, CS%Kd)
-  call uvchksum("[uv]h after ALE regridding/remapping of inputs", CS%uhtr, CS%vhtr, CS%G%HI)
+
+  if (CS%id_temp_regrid>0) call post_data(CS%id_temp_regrid, CS%tv%T, CS%diag, alt_h = CS%h_start)
+  if (CS%id_salt_regrid>0) call post_data(CS%id_salt_regrid, CS%tv%S, CS%diag, alt_h =CS%h_start)
+  if (CS%id_uhtr_regrid>0) call post_data(CS%id_uhtr_regrid, CS%uhtr, CS%diag, alt_h = CS%h_start)
+  if (CS%id_vhtr_regrid>0) call post_data(CS%id_vhtr_regrid, CS%vhtr, CS%diag, alt_h =CS%h_start)
+  if (CS%id_h_regrid>0) call post_data(CS%id_salt_regrid, h, CS%diag, alt_h = CS%h_start)
+
+  if (CS%debug) then
+    call uvchksum("[uv]h after ALE regridding/remapping of inputs", CS%uhtr, CS%vhtr, CS%G%HI)
+    call hchksum(CS%h_start,"h_start after update offline from files and arrays", CS%G%HI)
+  endif
 
   ! Update halos for some
   call pass_var(CS%h_end, CS%G%Domain) ! ALE remapping/regridding needs this
@@ -1061,6 +1090,19 @@ subroutine register_diags_offline_transport(Time, diag, CS)
   CS%id_h_redist = register_diag_field('ocean_model','h_redist', diag%axesTL, Time, &
     'Layer thicknesses before redistribution of mass fluxes','m')
 
+  ! Regridded/remapped input fields
+  CS%id_uhtr_regrid = register_diag_field('ocean_model', 'uhtr_regrid', diag%axesCuL, Time, &
+                                          'Zonal mass transport regridded/remapped onto offline grid','kg')
+  CS%id_vhtr_regrid = register_diag_field('ocean_model', 'vhtr_regrid', diag%axesCvL, Time, &
+                                          'Meridional mass transport regridded/remapped onto offline grid','kg')
+  CS%id_temp_regrid = register_diag_field('ocean_model', 'temp_regrid', diag%axesTL, Time, &
+                                          'Temperature regridded/remapped onto offline grid','C')
+  CS%id_salt_regrid = register_diag_field('ocean_model', 'salt_regrid', diag%axesTL, Time, &
+                                          'Salinity regridded/remapped onto offline grid','g kg-1')
+  CS%id_h_regrid = register_diag_field('ocean_model', 'h_regrid', diag%axesTL, Time, &
+                                          'Layer thicknesses regridded/remapped onto offline grid','m')
+
+
 end subroutine register_diags_offline_transport
 
 !> Posts diagnostics related to offline convergence diagnostics