diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index ad9e235b27..4b16730fee 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1180,7 +1180,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call cpu_clock_begin(id_clock_thermo) if (.not.CS%adiabatic) then if (CS%debug) then - call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2) + call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2, scale=US%L_T_to_m_s) call hchksum(h,"Pre-diabatic h", G%HI, haloshift=1, scale=GV%H_to_m) call uvchksum("Pre-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, & haloshift=0, scale=GV%H_to_m*US%L_to_m**2) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index c479550847..8c016b11b0 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -480,7 +480,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call disable_averaging(CS%diag) if (CS%debug) then - call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym) + call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s) endif call vertvisc_coef(up, vp, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp) @@ -670,7 +670,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (CS%debug) then call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym) - call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym) + call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_m) ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US) call check_redundant("Predictor up ", up, vp, G) @@ -860,7 +860,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag) if (CS%debug) then - call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym, vel_scale=1.0) + call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, scale=GV%H_to_m) ! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 05f2cac00a..9794070f20 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1284,20 +1284,22 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, ! surface mass flux maps handles%id_prcme = register_diag_field('ocean_model', 'PRCmE', diag%axesT1, Time, & - 'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1',& + 'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1', & standard_name='water_flux_into_sea_water', cmor_field_name='wfo', & cmor_standard_name='water_flux_into_sea_water',cmor_long_name='Water Flux Into Sea Water') + ! This diagnostic is rescaled to MKS units when combined. - handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time, & - 'Evaporation/condensation at ocean surface (evaporation is negative)', 'kg m-2 s-1',& - standard_name='water_evaporation_flux', cmor_field_name='evs', & - cmor_standard_name='water_evaporation_flux', & + handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time, & + 'Evaporation/condensation at ocean surface (evaporation is negative)', & + 'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & + standard_name='water_evaporation_flux', cmor_field_name='evs', & + cmor_standard_name='water_evaporation_flux', & cmor_long_name='Water Evaporation Flux Where Ice Free Ocean over Sea') ! smg: seaice_melt field requires updates to the sea ice model handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt', & diag%axesT1, Time, 'water flux to ocean from snow/sea ice melting(> 0) or formation(< 0)', & - 'kg m-2 s-1', & + 'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics', & cmor_field_name='fsitherm', & cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',& @@ -1305,6 +1307,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_precip = register_diag_field('ocean_model', 'precip', diag%axesT1, Time, & 'Liquid + frozen precipitation into ocean', 'kg m-2 s-1') + ! This diagnostic is rescaled to MKS units when combined. handles%id_fprec = register_diag_field('ocean_model', 'fprec', diag%axesT1, Time, & 'Frozen precipitation into ocean', & @@ -1324,32 +1327,39 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T) handles%id_frunoff = register_diag_field('ocean_model', 'frunoff', diag%axesT1, Time, & - 'Frozen runoff (calving) and iceberg melt into ocean', 'kg m-2 s-1', & + 'Frozen runoff (calving) and iceberg melt into ocean', & + units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & standard_name='water_flux_into_sea_water_from_icebergs', & cmor_field_name='ficeberg', & cmor_standard_name='water_flux_into_sea_water_from_icebergs', & cmor_long_name='Water Flux into Seawater from Icebergs') handles%id_lrunoff = register_diag_field('ocean_model', 'lrunoff', diag%axesT1, Time, & - 'Liquid runoff (rivers) into ocean', 'kg m-2 s-1', & + 'Liquid runoff (rivers) into ocean', & + units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & standard_name='water_flux_into_sea_water_from_rivers', cmor_field_name='friver', & cmor_standard_name='water_flux_into_sea_water_from_rivers', & cmor_long_name='Water Flux into Sea Water From Rivers') handles%id_net_massout = register_diag_field('ocean_model', 'net_massout', diag%axesT1, Time, & 'Net mass leaving the ocean due to evaporation, seaice formation', 'kg m-2 s-1') + ! This diagnostic is rescaled to MKS units when combined. handles%id_net_massin = register_diag_field('ocean_model', 'net_massin', diag%axesT1, Time, & 'Net mass entering ocean due to precip, runoff, ice melt', 'kg m-2 s-1') + ! This diagnostic is rescaled to MKS units when combined. handles%id_massout_flux = register_diag_field('ocean_model', 'massout_flux', diag%axesT1, Time, & 'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation)', & 'kg m-2', conversion=diag%GV%H_to_kg_m2) + ! This diagnostic is calculated in MKS units. handles%id_massin_flux = register_diag_field('ocean_model', 'massin_flux', diag%axesT1, Time, & 'Net mass flux of freshwater into the ocean (used in boundary flux calculation)', 'kg m-2') + ! This diagnostic is calculated in MKS units. + !========================================================================= - ! area integrated surface mass transport + ! area integrated surface mass transport, all are rescaled to MKS units before area integration. handles%id_total_prcme = register_scalar_field('ocean_model', 'total_PRCmE', Time, diag, & long_name='Area integrated net surface water flux (precip+melt+liq runoff+ice calving-evap)',& diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index f99b6d7f7c..668f185152 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -534,9 +534,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ call find_eta(h, tv, G, GV, US, eta) do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = (eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) + tmp1(i,j,k) = US%Z_to_m*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) enddo ; enddo ; enddo - vol_tot = US%Z_to_m*reproducing_sum(tmp1, sums=vol_lay) + vol_tot = reproducing_sum(tmp1, sums=vol_lay) + do k=1,nz ; vol_lay(k) = US%m_to_Z * vol_lay(k) ; enddo else do k=1,nz ; do j=js,je ; do i=is,ie tmp1(i,j,k) = H_to_kg_m2 * h(i,j,k) * areaTm(i,j) diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90 index 829afb851f..85d5ce452b 100644 --- a/src/framework/MOM_spatial_means.F90 +++ b/src/framework/MOM_spatial_means.F90 @@ -43,7 +43,7 @@ function global_area_mean(var, G, scale) do j=js,je ; do i=is,ie tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j)) enddo ; enddo - global_area_mean = reproducing_sum(tmpForSumming) * (G%US%m_to_L**2 * G%IareaT_global) + global_area_mean = reproducing_sum(tmpForSumming) * G%IareaT_global end function global_area_mean @@ -182,17 +182,20 @@ end function global_mass_integral !> Determine the global mean of a field along rows of constant i, returning it !! in a 1-d array using the local indexing. This uses reproducing sums. -subroutine global_i_mean(array, i_mean, G, mask, scale) +subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged real, dimension(SZJ_(G)), intent(out) :: i_mean !< Global mean of array along its i-axis real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: mask !< An array used for weighting the i-mean - real, optional, intent(in) :: scale !< A rescaling factor for the variable + optional, intent(in) :: mask !< An array used for weighting the i-mean + real, optional, intent(in) :: scale !< A rescaling factor for the output variable + real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal + !! calculations that is removed from the output ! Local variables type(EFP_type), allocatable, dimension(:) :: asum, mask_sum real :: scalefac ! A scaling factor for the variable. + real :: unscale ! A factor for undoing any internal rescaling before output. real :: mask_sum_r integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j @@ -201,6 +204,10 @@ subroutine global_i_mean(array, i_mean, G, mask, scale) idg_off = G%idg_offset ; jdg_off = G%jdg_offset scalefac = 1.0 ; if (present(scale)) scalefac = scale + unscale = 1.0 + if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then + scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale + endif ; endif call reset_EFP_overflow_error() allocate(asum(G%jsg:G%jeg)) @@ -253,24 +260,29 @@ subroutine global_i_mean(array, i_mean, G, mask, scale) enddo endif + if (unscale /= 1.0) then ; do j=js,je ; i_mean(j) = unscale*i_mean(j) ; enddo ; endif + deallocate(asum) end subroutine global_i_mean !> Determine the global mean of a field along rows of constant j, returning it !! in a 1-d array using the local indexing. This uses reproducing sums. -subroutine global_j_mean(array, j_mean, G, mask, scale) +subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged real, dimension(SZI_(G)), intent(out) :: j_mean !< Global mean of array along its j-axis real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: mask !< An array used for weighting the j-mean - real, optional, intent(in) :: scale !< A rescaling factor for the variable + optional, intent(in) :: mask !< An array used for weighting the j-mean + real, optional, intent(in) :: scale !< A rescaling factor for the output variable + real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal + !! calculations that is removed from the output ! Local variables type(EFP_type), allocatable, dimension(:) :: asum, mask_sum real :: mask_sum_r real :: scalefac ! A scaling factor for the variable. + real :: unscale ! A factor for undoing any internal rescaling before output. integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j @@ -278,6 +290,10 @@ subroutine global_j_mean(array, j_mean, G, mask, scale) idg_off = G%idg_offset ; jdg_off = G%jdg_offset scalefac = 1.0 ; if (present(scale)) scalefac = scale + unscale = 1.0 + if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then + scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale + endif ; endif call reset_EFP_overflow_error() allocate(asum(G%isg:G%ieg)) @@ -330,6 +346,8 @@ subroutine global_j_mean(array, j_mean, G, mask, scale) enddo endif + if (unscale /= 1.0) then ; do i=is,ie ; j_mean(i) = unscale*j_mean(i) ; enddo ; endif + deallocate(asum) end subroutine global_j_mean diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 8ed9a0a4c7..0ddca45c51 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -159,7 +159,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) call initialize_grid_rotation_angle(G, PF) ! Compute global integrals of grid values for later use in scalar diagnostics ! - call compute_global_grid_integrals(G) + call compute_global_grid_integrals(G, US=US) ! Write out all of the grid data used by this run. if (write_geom) call write_ocean_geometry_file(G, PF, output_dir, US=US) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 3d0fe6f1ed..3338f1fedb 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1145,17 +1145,21 @@ end subroutine set_velocity_depth_min ! ----------------------------------------------------------------------------- !> Pre-compute global integrals of grid quantities (like masked ocean area) for !! later use in reporting diagnostics -subroutine compute_global_grid_integrals(G) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid +subroutine compute_global_grid_integrals(G, US) + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming + real :: area_scale ! A scaling factor for area into MKS units integer :: i,j + area_scale = 1.0 ; if (present(US)) area_scale = US%L_to_m**2 + tmpForSumming(:,:) = 0. G%areaT_global = 0.0 ; G%IareaT_global = 0.0 do j=G%jsc,G%jec ; do i=G%isc,G%iec - tmpForSumming(i,j) = G%areaT(i,j) * G%mask2dT(i,j) + tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j) enddo ; enddo G%areaT_global = reproducing_sum(tmpForSumming) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 03310d70f3..ff08912191 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1889,16 +1889,19 @@ end subroutine set_velocity_depth_max !> Subroutine to pre-compute global integrals of grid quantities for !! later use in reporting diagnostics -subroutine compute_global_grid_integrals(G) +subroutine compute_global_grid_integrals(G, US) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming + real :: area_scale integer :: i,j + area_scale = US%L_to_m**2 tmpForSumming(:,:) = 0. G%areaT_global = 0.0 ; G%IareaT_global = 0.0 do j=G%jsc,G%jec ; do i=G%isc,G%iec - tmpForSumming(i,j) = G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) + tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j) enddo ; enddo G%areaT_global = reproducing_sum(tmpForSumming) G%IareaT_global = 1. / (G%areaT_global) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 38bf24ee60..a2257369a8 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -188,7 +188,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, scale=US%R_to_kg_m3*US%Z_to_m*US%L_to_m**2*US%s_to_T**3) if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, scale=US%L_T_to_m_s**2) call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, scale=US%s_to_T) - call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m) + call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m*US%L_to_m**2) endif sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index b4c100dc5d..eb1afb6bb8 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -3,6 +3,7 @@ module MOM_set_diffusivity ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_checksums, only : hchksum_pair use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type @@ -344,8 +345,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%useKappaShear) then if (CS%debug) then - call hchksum(u_h, "before calc_KS u_h",G%HI) - call hchksum(v_h, "before calc_KS v_h",G%HI) + call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, G%HI, scale=US%L_T_to_m_s) endif call cpu_clock_begin(id_clock_kappaShear) if (CS%Vertex_shear) then diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90 index dd0887845c..6016dbb98b 100644 --- a/src/parameterizations/vertical/MOM_sponge.F90 +++ b/src/parameterizations/vertical/MOM_sponge.F90 @@ -420,7 +420,7 @@ subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml) eta_anom(i,j) = e_D(i,j,k) - CS%Ref_eta_im(j,k) if (CS%Ref_eta_im(j,K) < -G%bathyT(i,j)) eta_anom(i,j) = 0.0 enddo ; enddo - call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G) + call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G, tmp_scale=US%Z_to_m) enddo if (CS%fldno > 0) allocate(fld_mean_anom(G%isd:G%ied,nz,CS%fldno)) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index e050933dc2..e425629c77 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -485,8 +485,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 + uhh(I)/(hprev(i+1,j,k)+h_neglect)) - CFL(I) = - uhh(I)/(hprev(i+1,j,k)+h_neglect) ! CFL is positive + !ts2(I) = 0.5*(1.0 + uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j))) + CFL(I) = - uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h hlos = MAX(0.0,-uhr(I-1,j,k)) @@ -497,8 +497,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect)) - CFL(I) = uhh(I)/(hprev(i,j,k)+h_neglect) ! CFL is positive + !ts2(I) = 0.5*(1.0 - uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) + CFL(I) = uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive endif enddo @@ -573,7 +573,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! Original implementation of PLM !flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I)) endif - !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect)) + !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j))) enddo ; enddo endif ! usePPM @@ -856,8 +856,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k)+h_neglect)) - CFL(i) = - vhh(i,J) / (hprev(i,j+1,k)+h_neglect) ! CFL is positive + !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) + CFL(i) = - vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h hlos = MAX(0.0,-vhr(i,J-1,k)) @@ -868,8 +868,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k)+h_neglect)) - CFL(i) = vhh(i,J) / (hprev(i,j,k)+h_neglect) ! CFL is positive + !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) + CFL(i) = vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive endif enddo