From 112ac4998c806077b3b1d2f7d3eab54d322347a5 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 11 Dec 2021 06:07:47 -0500 Subject: [PATCH] +(*)Revised offline tracer algorithms Minorly revised the algorithms used for offline tracer advection for rotational consistency, and for exact reproducibility across PE layouts by using reproducing sums to detect convergence. This also includes some changes to use roundoff to detect convergence instead of fixed values. Also replaced some divisions with multiplication by a reciprocal. In addition, some of the optional arguments to advect_tracer that are only used for offline tracer advection were renamed or revised for clarity and reordered for the convenience of the non-offline-tracer code. Although answers with the offline tracer code will change slightly because of this refactoring, because of some bugs that were fixed in another recent commit, it was previously impossible to have run the offline tracer cases at all. All answers and output in the MOM6-examples regression suite are bitwise identical. --- src/tracer/MOM_offline_aux.F90 | 132 +++++++++++++++---------------- src/tracer/MOM_offline_main.F90 | 129 +++++++++++++++--------------- src/tracer/MOM_tracer_advect.F90 | 49 ++++++------ 3 files changed, 154 insertions(+), 156 deletions(-) diff --git a/src/tracer/MOM_offline_aux.F90 b/src/tracer/MOM_offline_aux.F90 index b370dd6bb4..f95f2cd40e 100644 --- a/src/tracer/MOM_offline_aux.F90 +++ b/src/tracer/MOM_offline_aux.F90 @@ -70,7 +70,7 @@ subroutine update_h_horizontal_flux(G, GV, uhtr, vhtr, h_pre, h_new) max(GV%Angstrom_H, 1.0e-13*h_new(i,j,k) - G%areaT(i,j)*h_pre(i,j,k)) ! Convert back to thickness - h_new(i,j,k) = h_new(i,j,k) / (G%areaT(i,j)) + h_new(i,j,k) = h_new(i,j,k) * G%IareaT(i,j) enddo ; enddo enddo @@ -102,11 +102,11 @@ subroutine update_h_vertical_flux(G, GV, ea, eb, h_pre, h_new) do j=js-1,je+1 do i=is-1,ie+1 ! Top layer - h_new(i,j,1) = max(0.0, h_pre(i,j,1) + (eb(i,j,1) - ea(i,j,2) + ea(i,j,1) )) + h_new(i,j,1) = max(0.0, h_pre(i,j,1) + ((eb(i,j,1) - ea(i,j,2)) + ea(i,j,1))) h_new(i,j,1) = h_new(i,j,1) + max(0.0, 1.0e-13*h_new(i,j,1) - h_pre(i,j,1)) ! Bottom layer - h_new(i,j,nz) = max(0.0, h_pre(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1)+eb(i,j,nz))) + h_new(i,j,nz) = max(0.0, h_pre(i,j,nz) + ((ea(i,j,nz) - eb(i,j,nz-1)) + eb(i,j,nz))) h_new(i,j,nz) = h_new(i,j,nz) + max(0.0, 1.0e-13*h_new(i,j,nz) - h_pre(i,j,nz)) enddo @@ -140,13 +140,15 @@ subroutine limit_mass_flux_3d(G, GV, uh, vh, ea, eb, h_pre) !! step [H ~> m or kg m-2] ! Local variables - integer :: i, j, k, m, is, ie, js, je, nz - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: top_flux ! Net fluxes through the layer top [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: bottom_flux ! Net fluxes through the layer bottom [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: top_flux ! Net upward fluxes through the layer + ! top [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: bottom_flux ! Net downward fluxes through the layer + ! bottom [H ~> m or kg m-2] real :: pos_flux ! Net flux out of cell [H L2 ~> m3 or kg m-2] real :: hvol ! Cell volume [H L2 ~> m3 or kg m-2] real :: scale_factor ! A nondimensional rescaling factor between 0 and 1 [nondim] real :: max_off_cfl ! The maximum permitted fraction that can leave in a timestep [nondim] + integer :: i, j, k, m, is, ie, js, je, nz max_off_cfl = 0.5 @@ -182,46 +184,33 @@ subroutine limit_mass_flux_3d(G, GV, uh, vh, ea, eb, h_pre) do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 hvol = h_pre(i,j,k) * G%areaT(i,j) - pos_flux = max(0.0,-uh(I-1,j,k)) + max(0.0, -vh(i,J-1,k)) + & - max(0.0, uh(I,j,k)) + max(0.0, vh(i,J,k)) + & - max(0.0, top_flux(i,j,k)*G%areaT(i,j)) + max(0.0, bottom_flux(i,j,k)*G%areaT(i,j)) + pos_flux = ((max(0.0, -uh(I-1,j,k)) + max(0.0, uh(I,j,k))) + & + (max(0.0, -vh(i,J-1,k)) + max(0.0, vh(i,J,k)))) + & + (max(0.0, top_flux(i,j,k)) + max(0.0, bottom_flux(i,j,k))) * G%areaT(i,j) if (pos_flux>hvol .and. pos_flux>0.0) then - scale_factor = (hvol / pos_flux)*max_off_cfl + scale_factor = (hvol / pos_flux) * max_off_cfl else ! Don't scale scale_factor = 1.0 endif ! Scale horizontal fluxes - if (-uh(I-1,j,k)>0) uh(I-1,j,k) = uh(I-1,j,k)*scale_factor - if (uh(I,j,k)>0) uh(I,j,k) = uh(I,j,k)*scale_factor - if (-vh(i,J-1,k)>0) vh(i,J-1,k) = vh(i,J-1,k)*scale_factor - if (vh(i,J,k)>0) vh(i,J,k) = vh(i,J,k)*scale_factor - - if (k>1 .and. k0.0) then - ea(i,j,k) = ea(i,j,k)*scale_factor - eb(i,j,k-1) = eb(i,j,k-1)*scale_factor - endif - if (bottom_flux(i,j,k)>0.0) then - eb(i,j,k) = eb(i,j,k)*scale_factor - ea(i,j,k+1) = ea(i,j,k+1)*scale_factor - endif - ! Scale top layer - elseif (k==1) then - if (top_flux(i,j,k)>0.0) ea(i,j,k) = ea(i,j,k)*scale_factor - if (bottom_flux(i,j,k)>0.0) then - eb(i,j,k) = eb(i,j,k)*scale_factor - ea(i,j,k+1) = ea(i,j,k+1)*scale_factor - endif - ! Scale bottom layer - elseif (k==nz) then - if (top_flux(i,j,k)>0.0) then - ea(i,j,k) = ea(i,j,k)*scale_factor - eb(i,j,k-1) = eb(i,j,k-1)*scale_factor - endif - if (bottom_flux(i,j,k)>0.0) eb(i,j,k) = eb(i,j,k)*scale_factor + if (-uh(I-1,j,k) > 0.0) uh(I-1,j,k) = uh(I-1,j,k) * scale_factor + if (uh(I,j,k) > 0.0) uh(I,j,k) = uh(I,j,k) * scale_factor + if (-vh(i,J-1,k) > 0.0) vh(i,J-1,k) = vh(i,J-1,k) * scale_factor + if (vh(i,J,k) > 0.0) vh(i,J,k) = vh(i,J,k) * scale_factor + + ! Scale the flux across the interface atop a layer if it is upward + if (top_flux(i,j,k) > 0.0) then + ea(i,j,k) = ea(i,j,k) * scale_factor + if (k > 1) & + eb(i,j,k-1) = eb(i,j,k-1) * scale_factor + endif + ! Scale the flux across the interface atop a layer if it is downward + if (bottom_flux(i,j,k) > 0.0) then + eb(i,j,k) = eb(i,j,k) * scale_factor + if (k < nz) & + ea(i,j,k+1) = ea(i,j,k+1) * scale_factor endif enddo ; enddo ; enddo @@ -244,6 +233,8 @@ subroutine distribute_residual_uh_barotropic(G, GV, hvol, uh) real, dimension(SZI_(G),SZK_(GV)) :: h2d ! A 2-d slice of cell volumes [H L2 ~> m3 or kg] real, dimension(SZI_(G)) :: h2d_sum ! Vertically summed cell volumes [H L2 ~> m3 or kg] + real :: abs_uh_sum ! The vertical sum of the absolute value of the transports [H L2 ~> m3 or kg] + real :: new_uh_sum ! The vertically summed transports after redistribution [H L2 ~> m3 or kg] real :: uh_neglect ! A negligible transport [H L2 ~> m3 or kg] integer :: i, j, k, m, is, ie, js, je, nz @@ -253,7 +244,7 @@ subroutine distribute_residual_uh_barotropic(G, GV, hvol, uh) do j=js,je uh2d_sum(:) = 0.0 ! Copy over uh to a working array and sum up the remaining fluxes in a column - do k=1,nz ; do i=is-1,ie + do k=1,nz ; do I=is-1,ie uh2d(I,k) = uh(I,j,k) uh2d_sum(I) = uh2d_sum(I) + uh2d(I,k) enddo ; enddo @@ -265,13 +256,13 @@ subroutine distribute_residual_uh_barotropic(G, GV, hvol, uh) if (hvol(i,j,k)>0.) then h2d_sum(i) = h2d_sum(i) + h2d(i,k) else - h2d(i,k) = GV%H_subroundoff * 1.0*G%US%m_to_L**2 !### Change to G%areaT(i,j) + h2d(i,k) = GV%H_subroundoff * G%areaT(i,j) endif enddo ; enddo ! Distribute flux. Note min/max is intended to make sure that the mass transport ! does not deplete a cell - do i=is-1,ie + do I=is-1,ie if ( uh2d_sum(I)>0.0 ) then do k=1,nz uh2d(I,k) = uh2d_sum(I)*(h2d(i,k)/h2d_sum(i)) @@ -285,16 +276,20 @@ subroutine distribute_residual_uh_barotropic(G, GV, hvol, uh) uh2d(I,k) = 0.0 enddo endif - ! Calculate and check that column integrated transports match the original to - ! within the tolerance limit + + ! Check that column integrated transports match the original to within roundoff. uh_neglect = GV%Angstrom_H * min(G%areaT(i,j), G%areaT(i+1,j)) - ! ### This test may not work if GV%Angstrom_H is set to 0. - ! Instead try the max of this and ~roundoff compared with absolute transports? - if ( abs(sum(uh2d(I,:))-uh2d_sum(I)) > uh_neglect) & - call MOM_error(WARNING, "Column integral of uh does not match after barotropic redistribution") + abs_uh_sum = 0.0 ; new_uh_sum = 0.0 + do k=1,nz + abs_uh_sum = abs_uh_sum + abs(uh2d(j,k)) + new_uh_sum = new_uh_sum + uh2d(j,k) + enddo + if ( abs(new_uh_sum - uh2d_sum(j)) > max(uh_neglect, (5.0e-16*nz)*abs_uh_sum) ) & + call MOM_error(WARNING, "Column integral of uh does not match after "//& + "barotropic redistribution") enddo - do k=1,nz ; do i=is-1,ie + do k=1,nz ; do I=is-1,ie uh(I,j,k) = uh2d(I,k) enddo ; enddo enddo @@ -317,6 +312,8 @@ subroutine distribute_residual_vh_barotropic(G, GV, hvol, vh) real, dimension(SZJ_(G),SZK_(GV)) :: h2d ! A 2-d slice of cell volumes [H L2 ~> m3 or kg] real, dimension(SZJ_(G)) :: h2d_sum ! Vertically summed cell volumes [H L2 ~> m3 or kg] + real :: abs_vh_sum ! The vertical sum of the absolute value of the transports [H L2 ~> m3 or kg] + real :: new_vh_sum ! The vertically summed transports after redistribution [H L2 ~> m3 or kg] real :: vh_neglect ! A negligible transport [H L2 ~> m3 or kg] integer :: i, j, k, m, is, ie, js, je, nz @@ -326,7 +323,7 @@ subroutine distribute_residual_vh_barotropic(G, GV, hvol, vh) do i=is,ie vh2d_sum(:) = 0.0 ! Copy over uh to a working array and sum up the remaining fluxes in a column - do k=1,nz ; do j=js-1,je + do k=1,nz ; do J=js-1,je vh2d(J,k) = vh(i,J,k) vh2d_sum(J) = vh2d_sum(J) + vh2d(J,k) enddo ; enddo @@ -338,12 +335,12 @@ subroutine distribute_residual_vh_barotropic(G, GV, hvol, vh) if (hvol(i,j,k)>0.) then h2d_sum(j) = h2d_sum(j) + h2d(j,k) else - h2d(j,k) = GV%H_subroundoff * 1.0*G%US%m_to_L**2 !### Change to G%areaT(i,j) + h2d(i,k) = GV%H_subroundoff * G%areaT(i,j) endif enddo ; enddo ! Distribute flux evenly throughout a column - do j=js-1,je + do J=js-1,je if ( vh2d_sum(J)>0.0 ) then do k=1,nz vh2d(J,k) = vh2d_sum(J)*(h2d(j,k)/h2d_sum(j)) @@ -357,19 +354,20 @@ subroutine distribute_residual_vh_barotropic(G, GV, hvol, vh) vh2d(J,k) = 0.0 enddo endif - ! Calculate and check that column integrated transports match the original to - ! within the tolerance limit - vh_neglect = GV%Angstrom_H * min(G%areaT(i,j), G%areaT(i,j+1)) - ! ### This test may not work if GV%Angstrom_H is set to 0. - ! Instead try the max of this and ~roundoff compared with absolute transports? - if ( abs(sum(vh2d(J,:))-vh2d_sum(J)) > vh_neglect) then - call MOM_error(WARNING,"Column integral of vh does not match after "//& - "barotropic redistribution") - endif + ! Check that column integrated transports match the original to within roundoff. + vh_neglect = GV%Angstrom_H * min(G%areaT(i,j), G%areaT(i,j+1)) + abs_vh_sum = 0.0 ; new_vh_sum = 0.0 + do k=1,nz + abs_vh_sum = abs_vh_sum + abs(vh2d(J,k)) + new_vh_sum = new_vh_sum + vh2d(J,k) + enddo + if ( abs(new_vh_sum - vh2d_sum(J)) > max(vh_neglect, (5.0e-16*nz)*abs_vh_sum) ) & + call MOM_error(WARNING, "Column integral of vh does not match after "//& + "barotropic redistribution") enddo - do k=1,nz ; do j=js-1,je + do k=1,nz ; do J=js-1,je vh(i,J,k) = vh2d(J,k) enddo ; enddo enddo @@ -411,7 +409,7 @@ subroutine distribute_residual_uh_upwards(G, GV, hvol, uh) h2d(i,k) = hvol(i,j,k) - min_h * G%areaT(i,j) enddo ; enddo - do i=is-1,ie + do I=is-1,ie uh_col = SUM(uh2d(I,:)) ! Store original column-integrated transport do k=1,nz uh_remain = uh2d(I,k) @@ -466,7 +464,7 @@ subroutine distribute_residual_uh_upwards(G, GV, hvol, uh) enddo ! i-loop - do k=1,nz ; do i=is-1,ie + do k=1,nz ; do I=is-1,ie uh(I,j,k) = uh2d(I,k) enddo ; enddo enddo @@ -502,14 +500,14 @@ subroutine distribute_residual_vh_upwards(G, GV, hvol, vh) do i=is,ie ! Copy over uh and cell volume to working arrays - do k=1,nz ; do j=js-2,je+1 + do k=1,nz ; do J=js-2,je+1 vh2d(J,k) = vh(i,J,k) enddo ; enddo do k=1,nz ; do j=js-1,je+1 h2d(j,k) = hvol(i,j,k) - min_h * G%areaT(i,j) enddo ; enddo - do j=js-1,je + do J=js-1,je vh_col = SUM(vh2d(J,:)) do k=1,nz vh_remain = vh2d(J,k) @@ -565,7 +563,7 @@ subroutine distribute_residual_vh_upwards(G, GV, hvol, vh) endif enddo - do k=1,nz ; do j=js-1,je + do k=1,nz ; do J=js-1,je vh(i,J,k) = vh2d(J,k) enddo ; enddo enddo diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 72041fbc86..d5b3f708a3 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -6,6 +6,7 @@ module MOM_offline_main use MOM_ALE, only : ALE_CS, ALE_main_offline, ALE_offline_inputs 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 use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE @@ -13,7 +14,7 @@ module MOM_offline_main use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member use MOM_diabatic_aux, only : tridiagTS use MOM_diag_mediator, only : diag_ctrl, post_data, register_diag_field -use MOM_domains, only : sum_across_PEs, pass_var, pass_vector +use MOM_domains, only : pass_var, pass_vector use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : read_param, get_param, log_version, param_file_type @@ -39,7 +40,6 @@ module MOM_offline_main implicit none ; private #include "MOM_memory.h" -#include "version_variable.h" !> The control structure for the offline transport module type, public :: offline_transport_CS ; private @@ -305,7 +305,7 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C call hchksum(h_pre, "h_pre before transport", G%HI, scale=GV%H_to_m) call uvchksum("[uv]htr_sub before transport", uhtr_sub, vhtr_sub, G%HI, scale=HL2_to_kg_scale) endif - tot_residual = remaining_transport_sum(G, GV, uhtr, vhtr) + tot_residual = remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new) if (CS%print_adv_offline) then write(mesg,'(A,ES24.16)') "Main advection starting transport: ", tot_residual*HL2_to_kg_scale call MOM_mesg(mesg) @@ -328,15 +328,15 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C endif call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, CS%dt_offline, G, GV, US, & - CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=1, & - uhr_out=uhtr, vhr_out=vhtr, h_out=h_new, x_first_in=x_before_y) + CS%tracer_adv_CSp, CS%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, & + max_iter_in=1, update_vol_prev=.true., uhr_out=uhtr, vhr_out=vhtr) ! Switch the direction every iteration x_before_y = .not. x_before_y ! Update the new layer thicknesses after one round of advection has happened do k=1,nz ; do j=js,je ; do i=is,ie - h_new(i,j,k) = h_new(i,j,k) / (G%areaT(i,j)) !### Replace with "* G%IareaT(i,j)" + h_new(i,j,k) = h_vol(i,j,k) * G%IareaT(i,j) enddo ; enddo ; enddo if (MODULO(iter,CS%off_ale_mod)==0) then @@ -367,7 +367,7 @@ subroutine offline_advection_ale(fluxes, Time_start, time_interval, G, GV, US, C ! Check for whether we've used up all the advection, or if we need to move on because ! advection has stalled - tot_residual = remaining_transport_sum(G, GV, uhtr, vhtr) + tot_residual = remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new) if (CS%print_adv_offline) then write(mesg,'(A,ES24.16)') "Main advection remaining transport: ", tot_residual*HL2_to_kg_scale call MOM_mesg(mesg) @@ -478,9 +478,6 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve call pass_var(h_vol,G%Domain) call pass_vector(uhtr, vhtr, G%Domain) - ! Store volumes for advect_tracer - h_pre(:,:,:) = h_vol(:,:,:) - if (CS%debug) then call MOM_tracer_chksum("Before upwards redistribute ", CS%tracer_Reg%Tr, CS%tracer_Reg%ntr, G) call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, G%HI, scale=HL2_to_kg_scale) @@ -495,8 +492,8 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve endif call advect_tracer(h_pre, uhtr, vhtr, CS%OBC, CS%dt_offline, G, GV, US, & - CS%tracer_adv_CSp, CS%tracer_Reg, h_prev_opt=h_pre, max_iter_in=1, & - h_out=h_vol, uhr_out=uhr, vhr_out=vhr, x_first_in=x_before_y) + CS%tracer_adv_CSp, CS%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, & + max_iter_in=1, update_vol_prev=.true., uhr_out=uhr, vhr_out=vhr) if (CS%debug) then call MOM_tracer_chksum("After upwards redistribute ", CS%tracer_Reg%Tr, CS%tracer_Reg%ntr, G) @@ -506,7 +503,7 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve do k=1,nz ; do j=js,je ; do i=is,ie uhtr(I,j,k) = uhr(I,j,k) vhtr(i,J,k) = vhr(i,J,k) - h_new(i,j,k) = h_vol(i,j,k) / (G%areaT(i,j)) + h_new(i,j,k) = h_vol(i,j,k) * G%IareaT(i,j) h_pre(i,j,k) = h_new(i,j,k) enddo ; enddo ; enddo @@ -522,9 +519,6 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve call pass_var(h_vol, G%Domain) call pass_vector(uhtr, vhtr, G%Domain) - ! Copy h_vol to h_pre for advect_tracer routine - h_pre(:,:,:) = h_vol(:,:,:) - if (CS%debug) then call MOM_tracer_chksum("Before barotropic redistribute ", CS%tracer_Reg%Tr, CS%tracer_Reg%ntr, G) call uvchksum("[uv]tr before upwards redistribute", uhtr, vhtr, G%HI, scale=HL2_to_kg_scale) @@ -539,8 +533,8 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve endif call advect_tracer(h_pre, uhtr, vhtr, CS%OBC, CS%dt_offline, G, GV, US, & - CS%tracer_adv_CSp, CS%tracer_Reg, h_prev_opt=h_pre, max_iter_in=1, & - h_out=h_vol, uhr_out=uhr, vhr_out=vhr, x_first_in=x_before_y) + CS%tracer_adv_CSp, CS%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, & + max_iter_in=1, update_vol_prev=.true., uhr_out=uhr, vhr_out=vhr) if (CS%debug) then call MOM_tracer_chksum("After barotropic redistribute ", CS%tracer_Reg%Tr, CS%tracer_Reg%ntr, G) @@ -550,14 +544,14 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve do k=1,nz ; do j=js,je ; do i=is,ie uhtr(I,j,k) = uhr(I,j,k) vhtr(i,J,k) = vhr(i,J,k) - h_new(i,j,k) = h_vol(i,j,k) / (G%areaT(i,j)) + h_new(i,j,k) = h_vol(i,j,k) * G%IareaT(i,j) h_pre(i,j,k) = h_new(i,j,k) enddo ; enddo ; enddo endif ! redistribute barotropic ! Check to see if all transport has been exhausted - tot_residual = remaining_transport_sum(G, GV, uhtr, vhtr) + tot_residual = remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new) if (CS%print_adv_offline) then write(mesg,'(A,ES24.16)') "Residual advection remaining transport: ", tot_residual*HL2_to_kg_scale call MOM_mesg(mesg) @@ -597,39 +591,40 @@ subroutine offline_redistribute_residual(CS, G, GV, US, h_pre, uhtr, vhtr, conve end subroutine offline_redistribute_residual !> Returns the sums of any non-negligible remaining transport [H L2 ~> m3 or kg] to check for advection convergence -real function remaining_transport_sum(G, GV, uhtr, vhtr) +real function remaining_transport_sum(G, GV, US, uhtr, vhtr, h_new) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(in ) :: uhtr !< Zonal mass transport [H L2 ~> m3 or kg] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(in ) :: vhtr !< Meridional mass transport [H L2 ~> m3 or kg] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in ) :: h_new !< Layer thicknesses [H ~> m or kg m-2] ! Local variables - integer :: i, j, k - integer :: is, ie, js, je, nz - real :: h_min !< A layer thickness below roundoff from GV type - real :: uh_neglect !< A small value of zonal transport that effectively is below roundoff error - real :: vh_neglect !< A small value of meridional transport that effectively is below roundoff error + real, dimension(SZI_(G),SZJ_(G)) :: trans_rem_col !< The vertical sum of the absolute value of + !! transports through the faces of a column, in MKS units [kg]. + real :: trans_cell !< The sum of the absolute value of the remaining transports through the faces + !! of a tracer cell [H L2 ~> m3 or kg] + real :: HL2_to_kg_scale !< Unit conversion factor to cell mass [kg H-1 L-2 ~> kg m-3 or 1] + integer :: i, j, k, is, ie, js, je, nz - nz = GV%ke - 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 - h_min = GV%H_subroundoff + HL2_to_kg_scale = GV%H_to_kg_m2 * US%L_to_m**2 - remaining_transport_sum = 0. + trans_rem_col(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie - uh_neglect = h_min * MIN(G%areaT(i,j), G%areaT(i+1,j)) - vh_neglect = h_min * MIN(G%areaT(i,j), G%areaT(i,j+1)) - if (ABS(uhtr(I,j,k))>uh_neglect) then - remaining_transport_sum = remaining_transport_sum + ABS(uhtr(I,j,k)) - endif - if (ABS(vhtr(i,J,k))>vh_neglect) then - remaining_transport_sum = remaining_transport_sum + ABS(vhtr(i,J,k)) - endif + trans_cell = (ABS(uhtr(I-1,j,k)) + ABS(uhtr(I,j,k))) + & + (ABS(vhtr(i,J-1,k)) + ABS(vhtr(i,J,k))) + if (trans_cell > max(1.0e-16*h_new(i,j,k), GV%H_subroundoff) * G%areaT(i,j)) & + trans_rem_col(i,j) = trans_rem_col(i,j) + HL2_to_kg_scale * trans_cell enddo ; enddo ; enddo - !### The value of this sum is not layout independent. - call sum_across_PEs(remaining_transport_sum) + + ! The factor of 0.5 here is to avoid double-counting because two cells share a face. + remaining_transport_sum = 0.5 * GV%kg_m2_to_H*US%m_to_L**2 * & + reproducing_sum(trans_rem_col, is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd)) end function remaining_transport_sum @@ -854,11 +849,14 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, ! Remaining meridional mass transports [H L2 ~> m3 or kg] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vhtr_sub - real :: sum_abs_fluxes, sum_u, sum_v ! Used to keep track of how close to convergence we are [H L2 ~> m3 or kg] - ! Vertical diffusion related variables [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJB_(G)) :: rem_col_flux ! The summed absolute value of the remaining + ! fluxes through the faces of a column or within a column, in mks units [kg] + real :: sum_flux ! Globally summed absolute value of fluxes in mks units [kg], which is + ! used to keep track of how close to convergence we are. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & - eatr_sub, & - ebtr_sub + eatr_sub, & ! Layer entrainment rate from above for this sub-cycle [H ~> m or kg m-2] + ebtr_sub ! Layer entrainment rate from below for this sub-cycle [H ~> m or kg m-2] ! 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 thicknesses [H ~> m or kg m-2] @@ -899,7 +897,6 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, vhtr_sub(i,J,k) = vhtr(i,J,k) enddo ; enddo ; enddo - ! Calculate 3d mass transports to be used in this iteration call limit_mass_flux_3d(G, GV, uhtr_sub, vhtr_sub, eatr_sub, ebtr_sub, h_pre) @@ -920,11 +917,11 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, h_vol(i,j,k) = h_pre(i,j,k) * G%areaT(i,j) enddo ; enddo ; enddo call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, US, & - CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=30, x_first_in=x_before_y) + CS%tracer_adv_CSp, CS%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, max_iter_in=30) ! Done with horizontal so now h_pre should be h_new do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1 - h_pre(i,j,k) = h_new(i,j,k) + h_pre(i,j,k) = h_new(i,j,k) enddo ; enddo ; enddo endif @@ -936,12 +933,12 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1 h_vol(i,j,k) = h_pre(i,j,k) * G%areaT(i,j) enddo ; enddo ; enddo - call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, US, & - CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=30, x_first_in=x_before_y) + call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, US, CS%tracer_adv_CSp, & + CS%tracer_Reg, x_first_in=x_before_y, vol_prev=h_vol, max_iter_in=30) ! Done with horizontal so now h_pre should be h_new do k=1,nz ; do i=is-1,ie+1 ; do j=js-1,je+1 - h_pre(i,j,k) = h_new(i,j,k) + h_pre(i,j,k) = h_new(i,j,k) enddo ; enddo ; enddo ! Second vertical advection @@ -973,28 +970,25 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, G, GV, US, call pass_var(ebtr,G%Domain) call pass_var(h_pre,G%Domain) call pass_vector(uhtr,vhtr,G%Domain) - ! + ! Calculate how close we are to converging by summing the remaining fluxes at each point - sum_abs_fluxes = 0.0 - sum_u = 0.0 - sum_v = 0.0 + HL2_to_kg_scale = US%L_to_m**2*GV%H_to_kg_m2 + rem_col_flux(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie - sum_u = sum_u + abs(uhtr(I-1,j,k))+abs(uhtr(I,j,k)) - sum_v = sum_v + abs(vhtr(i,J-1,k))+abs(vhtr(I,J,k)) - sum_abs_fluxes = sum_abs_fluxes + abs(eatr(i,j,k)) + abs(ebtr(i,j,k)) + abs(uhtr(I-1,j,k)) + & - abs(uhtr(I,j,k)) + abs(vhtr(i,J-1,k)) + abs(vhtr(i,J,k)) + rem_col_flux(i,j) = rem_col_flux(i,j) + HL2_to_kg_scale * & + ( (abs(eatr(i,j,k)) + abs(ebtr(i,j,k))) + & + ((abs(uhtr(I-1,j,k)) + abs(uhtr(I,j,k))) + & + (abs(vhtr(i,J-1,k)) + abs(vhtr(i,J,k))) ) ) enddo ; enddo ; enddo - call sum_across_PEs(sum_abs_fluxes) - - HL2_to_kg_scale = US%L_to_m**2*GV%H_to_kg_m2 + sum_flux = reproducing_sum(rem_col_flux, is+(1-G%isd), ie+(1-G%isd), js+(1-G%jsd), je+(1-G%jsd)) - write(mesg,*) "offline_advection_layer: Remaining u-flux, v-flux:", & - sum_u*HL2_to_kg_scale, sum_v*HL2_to_kg_scale - call MOM_mesg(mesg) - if (sum_abs_fluxes==0) then + if (sum_flux==0) then write(mesg,*) 'offline_advection_layer: Converged after iteration', iter call MOM_mesg(mesg) exit + else + write(mesg,*) "offline_advection_layer: Iteration ", iter, " remaining total fluxes: ", sum_flux + call MOM_mesg(mesg) endif ! Switch order of Strang split every iteration @@ -1321,7 +1315,8 @@ subroutine offline_transport_init(param_file, CS, diabatic_CSp, G, GV, US) character(len=40) :: mdl = "offline_transport" character(len=20) :: redistribute_method - + ! This include declares and sets the variable "version". +# include "version_variable.h" integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB @@ -1336,7 +1331,7 @@ subroutine offline_transport_init(param_file, CS, diabatic_CSp, G, GV, US) return endif allocate(CS) - call log_version(param_file, mdl,version, "This module allows for tracers to be run offline") + call log_version(param_file, mdl, version, "This module allows for tracers to be run offline") ! Parse MOM_input for offline control call get_param(param_file, mdl, "OFFLINEDIR", CS%offlinedir, & diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 1ad6343cf8..e2c669fcc7 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -47,36 +47,41 @@ module MOM_tracer_advect !> This routine time steps the tracer concentration using a !! monotonic, conservative, weakly diffusive scheme. -subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & - h_prev_opt, max_iter_in, x_first_in, uhr_out, vhr_out, h_out) +subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first_in, & + vol_prev, max_iter_in, update_vol_prev, uhr_out, vhr_out) type(ocean_grid_type), intent(inout) :: 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_end !< layer thickness after advection [H ~> m or kg m-2] + intent(in) :: h_end !< Layer thickness after advection [H ~> m or kg m-2] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: uhtr !< accumulated volume/mass flux through zonal face [H L2 ~> m3 or kg] + intent(in) :: uhtr !< Accumulated volume or mass flux through the + !! zonal faces [H L2 ~> m3 or kg] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(in) :: vhtr !< accumulated volume/mass flux through merid face [H L2 ~> m3 or kg] + intent(in) :: vhtr !< Accumulated volume or mass flux through the + !! meridional faces [H L2 ~> m3 or kg] type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used real, intent(in) :: dt !< time increment [T ~> s] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(tracer_advect_CS), pointer :: CS !< control structure for module type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - optional, intent(in) :: h_prev_opt !< Cell volume before advection [H L2 ~> m3 or kg] - integer, optional, intent(in) :: max_iter_in !< The maximum number of iterations logical, optional, intent(in) :: x_first_in !< If present, indicate whether to update !! first in the x- or y-direction. + ! The remaining optional arguments are only used in offline tracer mode. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + optional, intent(inout) :: vol_prev !< Cell volume before advection [H L2 ~> m3 or kg]. + !! If update_vol_prev is true, the returned value is + !! the cell volume after the transport that was done + !! by this call, and if all the transport could be + !! accommodated it should be close to h_end*G%areaT. + integer, optional, intent(in) :: max_iter_in !< The maximum number of iterations + logical, optional, intent(in) :: update_vol_prev !< If present and true, update vol_prev to + !! return its value after the tracer have been updated. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - optional, intent(out) :: uhr_out !< Remaining accumulated volume/mass flux through zonal face - !! [H L2 ~> m3 or kg] + optional, intent(out) :: uhr_out !< Remaining accumulated volume or mass fluxes + !! through the zonal faces [H L2 ~> m3 or kg] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - optional, intent(out) :: vhr_out !< Remaining accumulated volume/mass flux through meridional face - !! [H L2 ~> m3 or kg] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - optional, intent(out) :: h_out !< Cell volume after the transport that was done - !! by this call [H L2 ~> m3 or kg]. If all the transport - !! could be accommodated, this is close to h_end*G%areaT. + optional, intent(out) :: vhr_out !< Remaining accumulated volume or mass fluxes + !! through the meridional faces [H L2 ~> m3 or kg] type(tracer_type) :: Tr(MAX_FIELDS_) ! The array of registered tracers real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & @@ -137,9 +142,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & enddo call cpu_clock_end(id_clock_pass) -!$OMP parallel default(none) shared(nz,jsd,jed,IsdB,IedB,uhr,jsdB,jedB,Isd,Ied,vhr, & -!$OMP hprev,domore_k,js,je,is,ie,uhtr,vhtr,G,GV,h_end,& -!$OMP uh_neglect,vh_neglect,ntr,Tr,h_prev_opt) + !$OMP parallel default(shared) ! This initializes the halos of uhr and vhr because pass_vector might do ! calculations on them, even though they are never used. @@ -152,7 +155,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & ! Put the remaining (total) thickness fluxes into uhr and vhr. do j=js,je ; do I=is-1,ie ; uhr(I,j,k) = uhtr(I,j,k) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; vhr(i,J,k) = vhtr(i,J,k) ; enddo ; enddo - if (.not. present(h_prev_opt)) then + if (.not. present(vol_prev)) then ! This loop reconstructs the thickness field the last time that the ! tracers were updated, probably just after the diabatic forcing. A useful ! diagnostic could be to compare this reconstruction with that older value. @@ -167,7 +170,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & enddo ; enddo else do i=is,ie ; do j=js,je - hprev(i,j,k) = h_prev_opt(i,j,k) + hprev(i,j,k) = vol_prev(i,j,k) enddo ; enddo endif enddo @@ -326,7 +329,9 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & if (present(uhr_out)) uhr_out(:,:,:) = uhr(:,:,:) if (present(vhr_out)) vhr_out(:,:,:) = vhr(:,:,:) - if (present(h_out)) h_out(:,:,:) = hprev(:,:,:) + if (present(vol_prev) .and. present(update_vol_prev)) then + if (update_vol_prev) vol_prev(:,:,:) = hprev(:,:,:) + endif call cpu_clock_end(id_clock_advect)