Skip to content

Commit 4d7a947

Browse files
marshallwardmarshallward
authored
Merge pull request mom-ocean#1040 from Hallberg-NOAA/rescale_diag_fix
MOM6: +(*)Dimensional consistency completion
2 parents c9a7f98 + 274a613 commit 4d7a947

12 files changed

+77
-41
lines changed

‎src/core/MOM.F90

+1-1
Original file line numberDiff line numberDiff line change
@@ -1180,7 +1180,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
11801180
call cpu_clock_begin(id_clock_thermo)
11811181
if (.not.CS%adiabatic) then
11821182
if (CS%debug) then
1183-
call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2)
1183+
call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2, scale=US%L_T_to_m_s)
11841184
call hchksum(h,"Pre-diabatic h", G%HI, haloshift=1, scale=GV%H_to_m)
11851185
call uvchksum("Pre-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, &
11861186
haloshift=0, scale=GV%H_to_m*US%L_to_m**2)

‎src/core/MOM_dynamics_split_RK2.F90

+3-3
Original file line numberDiff line numberDiff line change
@@ -480,7 +480,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
480480
call disable_averaging(CS%diag)
481481

482482
if (CS%debug) then
483-
call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym)
483+
call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
484484
endif
485485
call vertvisc_coef(up, vp, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC)
486486
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, &
670670

671671
if (CS%debug) then
672672
call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym)
673-
call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym)
673+
call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
674674
call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_m)
675675
! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
676676
call check_redundant("Predictor up ", up, vp, G)
@@ -860,7 +860,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
860860
if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag)
861861

862862
if (CS%debug) then
863-
call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym, vel_scale=1.0)
863+
call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
864864
call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
865865
call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, scale=GV%H_to_m)
866866
! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)

‎src/core/MOM_forcing_type.F90

+19-9
Original file line numberDiff line numberDiff line change
@@ -1284,27 +1284,30 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
12841284
! surface mass flux maps
12851285

12861286
handles%id_prcme = register_diag_field('ocean_model', 'PRCmE', diag%axesT1, Time, &
1287-
'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1',&
1287+
'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1', &
12881288
standard_name='water_flux_into_sea_water', cmor_field_name='wfo', &
12891289
cmor_standard_name='water_flux_into_sea_water',cmor_long_name='Water Flux Into Sea Water')
1290+
! This diagnostic is rescaled to MKS units when combined.
12901291

1291-
handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time, &
1292-
'Evaporation/condensation at ocean surface (evaporation is negative)', 'kg m-2 s-1',&
1293-
standard_name='water_evaporation_flux', cmor_field_name='evs', &
1294-
cmor_standard_name='water_evaporation_flux', &
1292+
handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time, &
1293+
'Evaporation/condensation at ocean surface (evaporation is negative)', &
1294+
'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
1295+
standard_name='water_evaporation_flux', cmor_field_name='evs', &
1296+
cmor_standard_name='water_evaporation_flux', &
12951297
cmor_long_name='Water Evaporation Flux Where Ice Free Ocean over Sea')
12961298

12971299
! smg: seaice_melt field requires updates to the sea ice model
12981300
handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt', &
12991301
diag%axesT1, Time, 'water flux to ocean from snow/sea ice melting(> 0) or formation(< 0)', &
1300-
'kg m-2 s-1', &
1302+
'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
13011303
standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics', &
13021304
cmor_field_name='fsitherm', &
13031305
cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',&
13041306
cmor_long_name='water flux to ocean from sea ice melt(> 0) or form(< 0)')
13051307

13061308
handles%id_precip = register_diag_field('ocean_model', 'precip', diag%axesT1, Time, &
13071309
'Liquid + frozen precipitation into ocean', 'kg m-2 s-1')
1310+
! This diagnostic is rescaled to MKS units when combined.
13081311

13091312
handles%id_fprec = register_diag_field('ocean_model', 'fprec', diag%axesT1, Time, &
13101313
'Frozen precipitation into ocean', &
@@ -1324,32 +1327,39 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
13241327
units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T)
13251328

13261329
handles%id_frunoff = register_diag_field('ocean_model', 'frunoff', diag%axesT1, Time, &
1327-
'Frozen runoff (calving) and iceberg melt into ocean', 'kg m-2 s-1', &
1330+
'Frozen runoff (calving) and iceberg melt into ocean', &
1331+
units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
13281332
standard_name='water_flux_into_sea_water_from_icebergs', &
13291333
cmor_field_name='ficeberg', &
13301334
cmor_standard_name='water_flux_into_sea_water_from_icebergs', &
13311335
cmor_long_name='Water Flux into Seawater from Icebergs')
13321336

13331337
handles%id_lrunoff = register_diag_field('ocean_model', 'lrunoff', diag%axesT1, Time, &
1334-
'Liquid runoff (rivers) into ocean', 'kg m-2 s-1', &
1338+
'Liquid runoff (rivers) into ocean', &
1339+
units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
13351340
standard_name='water_flux_into_sea_water_from_rivers', cmor_field_name='friver', &
13361341
cmor_standard_name='water_flux_into_sea_water_from_rivers', &
13371342
cmor_long_name='Water Flux into Sea Water From Rivers')
13381343

13391344
handles%id_net_massout = register_diag_field('ocean_model', 'net_massout', diag%axesT1, Time, &
13401345
'Net mass leaving the ocean due to evaporation, seaice formation', 'kg m-2 s-1')
1346+
! This diagnostic is rescaled to MKS units when combined.
13411347

13421348
handles%id_net_massin = register_diag_field('ocean_model', 'net_massin', diag%axesT1, Time, &
13431349
'Net mass entering ocean due to precip, runoff, ice melt', 'kg m-2 s-1')
1350+
! This diagnostic is rescaled to MKS units when combined.
13441351

13451352
handles%id_massout_flux = register_diag_field('ocean_model', 'massout_flux', diag%axesT1, Time, &
13461353
'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation)', &
13471354
'kg m-2', conversion=diag%GV%H_to_kg_m2)
1355+
! This diagnostic is calculated in MKS units.
13481356

13491357
handles%id_massin_flux = register_diag_field('ocean_model', 'massin_flux', diag%axesT1, Time, &
13501358
'Net mass flux of freshwater into the ocean (used in boundary flux calculation)', 'kg m-2')
1359+
! This diagnostic is calculated in MKS units.
1360+
13511361
!=========================================================================
1352-
! area integrated surface mass transport
1362+
! area integrated surface mass transport, all are rescaled to MKS units before area integration.
13531363

13541364
handles%id_total_prcme = register_scalar_field('ocean_model', 'total_PRCmE', Time, diag, &
13551365
long_name='Area integrated net surface water flux (precip+melt+liq runoff+ice calving-evap)',&

‎src/diagnostics/MOM_sum_output.F90

+3-2
Original file line numberDiff line numberDiff line change
@@ -534,9 +534,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
534534

535535
call find_eta(h, tv, G, GV, US, eta)
536536
do k=1,nz ; do j=js,je ; do i=is,ie
537-
tmp1(i,j,k) = (eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j)
537+
tmp1(i,j,k) = US%Z_to_m*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j)
538538
enddo ; enddo ; enddo
539-
vol_tot = US%Z_to_m*reproducing_sum(tmp1, sums=vol_lay)
539+
vol_tot = reproducing_sum(tmp1, sums=vol_lay)
540+
do k=1,nz ; vol_lay(k) = US%m_to_Z * vol_lay(k) ; enddo
540541
else
541542
do k=1,nz ; do j=js,je ; do i=is,ie
542543
tmp1(i,j,k) = H_to_kg_m2 * h(i,j,k) * areaTm(i,j)

‎src/framework/MOM_spatial_means.F90

+25-7
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@ function global_area_mean(var, G, scale)
4343
do j=js,je ; do i=is,ie
4444
tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j))
4545
enddo ; enddo
46-
global_area_mean = reproducing_sum(tmpForSumming) * (G%US%m_to_L**2 * G%IareaT_global)
46+
global_area_mean = reproducing_sum(tmpForSumming) * G%IareaT_global
4747

4848
end function global_area_mean
4949

@@ -182,17 +182,20 @@ end function global_mass_integral
182182

183183
!> Determine the global mean of a field along rows of constant i, returning it
184184
!! in a 1-d array using the local indexing. This uses reproducing sums.
185-
subroutine global_i_mean(array, i_mean, G, mask, scale)
185+
subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale)
186186
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
187187
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged
188188
real, dimension(SZJ_(G)), intent(out) :: i_mean !< Global mean of array along its i-axis
189189
real, dimension(SZI_(G),SZJ_(G)), &
190-
optional, intent(in) :: mask !< An array used for weighting the i-mean
191-
real, optional, intent(in) :: scale !< A rescaling factor for the variable
190+
optional, intent(in) :: mask !< An array used for weighting the i-mean
191+
real, optional, intent(in) :: scale !< A rescaling factor for the output variable
192+
real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal
193+
!! calculations that is removed from the output
192194

193195
! Local variables
194196
type(EFP_type), allocatable, dimension(:) :: asum, mask_sum
195197
real :: scalefac ! A scaling factor for the variable.
198+
real :: unscale ! A factor for undoing any internal rescaling before output.
196199
real :: mask_sum_r
197200
integer :: is, ie, js, je, idg_off, jdg_off
198201
integer :: i, j
@@ -201,6 +204,10 @@ subroutine global_i_mean(array, i_mean, G, mask, scale)
201204
idg_off = G%idg_offset ; jdg_off = G%jdg_offset
202205

203206
scalefac = 1.0 ; if (present(scale)) scalefac = scale
207+
unscale = 1.0
208+
if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then
209+
scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale
210+
endif ; endif
204211
call reset_EFP_overflow_error()
205212

206213
allocate(asum(G%jsg:G%jeg))
@@ -253,31 +260,40 @@ subroutine global_i_mean(array, i_mean, G, mask, scale)
253260
enddo
254261
endif
255262

263+
if (unscale /= 1.0) then ; do j=js,je ; i_mean(j) = unscale*i_mean(j) ; enddo ; endif
264+
256265
deallocate(asum)
257266

258267
end subroutine global_i_mean
259268

260269
!> Determine the global mean of a field along rows of constant j, returning it
261270
!! in a 1-d array using the local indexing. This uses reproducing sums.
262-
subroutine global_j_mean(array, j_mean, G, mask, scale)
271+
subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale)
263272
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
264273
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged
265274
real, dimension(SZI_(G)), intent(out) :: j_mean !< Global mean of array along its j-axis
266275
real, dimension(SZI_(G),SZJ_(G)), &
267-
optional, intent(in) :: mask !< An array used for weighting the j-mean
268-
real, optional, intent(in) :: scale !< A rescaling factor for the variable
276+
optional, intent(in) :: mask !< An array used for weighting the j-mean
277+
real, optional, intent(in) :: scale !< A rescaling factor for the output variable
278+
real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal
279+
!! calculations that is removed from the output
269280

270281
! Local variables
271282
type(EFP_type), allocatable, dimension(:) :: asum, mask_sum
272283
real :: mask_sum_r
273284
real :: scalefac ! A scaling factor for the variable.
285+
real :: unscale ! A factor for undoing any internal rescaling before output.
274286
integer :: is, ie, js, je, idg_off, jdg_off
275287
integer :: i, j
276288

277289
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
278290
idg_off = G%idg_offset ; jdg_off = G%jdg_offset
279291

280292
scalefac = 1.0 ; if (present(scale)) scalefac = scale
293+
unscale = 1.0
294+
if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then
295+
scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale
296+
endif ; endif
281297
call reset_EFP_overflow_error()
282298

283299
allocate(asum(G%isg:G%ieg))
@@ -330,6 +346,8 @@ subroutine global_j_mean(array, j_mean, G, mask, scale)
330346
enddo
331347
endif
332348

349+
if (unscale /= 1.0) then ; do i=is,ie ; j_mean(i) = unscale*j_mean(i) ; enddo ; endif
350+
333351
deallocate(asum)
334352

335353
end subroutine global_j_mean

‎src/initialization/MOM_fixed_initialization.F90

+1-1
Original file line numberDiff line numberDiff line change
@@ -159,7 +159,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
159159
call initialize_grid_rotation_angle(G, PF)
160160

161161
! Compute global integrals of grid values for later use in scalar diagnostics !
162-
call compute_global_grid_integrals(G)
162+
call compute_global_grid_integrals(G, US=US)
163163

164164
! Write out all of the grid data used by this run.
165165
if (write_geom) call write_ocean_geometry_file(G, PF, output_dir, US=US)

‎src/initialization/MOM_shared_initialization.F90

+7-3
Original file line numberDiff line numberDiff line change
@@ -1145,17 +1145,21 @@ end subroutine set_velocity_depth_min
11451145
! -----------------------------------------------------------------------------
11461146
!> Pre-compute global integrals of grid quantities (like masked ocean area) for
11471147
!! later use in reporting diagnostics
1148-
subroutine compute_global_grid_integrals(G)
1149-
type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
1148+
subroutine compute_global_grid_integrals(G, US)
1149+
type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
1150+
type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
11501151

11511152
! Local variables
11521153
real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
1154+
real :: area_scale ! A scaling factor for area into MKS units
11531155
integer :: i,j
11541156

1157+
area_scale = 1.0 ; if (present(US)) area_scale = US%L_to_m**2
1158+
11551159
tmpForSumming(:,:) = 0.
11561160
G%areaT_global = 0.0 ; G%IareaT_global = 0.0
11571161
do j=G%jsc,G%jec ; do i=G%isc,G%iec
1158-
tmpForSumming(i,j) = G%areaT(i,j) * G%mask2dT(i,j)
1162+
tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j)
11591163
enddo ; enddo
11601164
G%areaT_global = reproducing_sum(tmpForSumming)
11611165

‎src/initialization/MOM_state_initialization.F90

+5-2
Original file line numberDiff line numberDiff line change
@@ -1889,16 +1889,19 @@ end subroutine set_velocity_depth_max
18891889

18901890
!> Subroutine to pre-compute global integrals of grid quantities for
18911891
!! later use in reporting diagnostics
1892-
subroutine compute_global_grid_integrals(G)
1892+
subroutine compute_global_grid_integrals(G, US)
18931893
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
1894+
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
18941895
! Local variables
18951896
real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
1897+
real :: area_scale
18961898
integer :: i,j
18971899

1900+
area_scale = US%L_to_m**2
18981901
tmpForSumming(:,:) = 0.
18991902
G%areaT_global = 0.0 ; G%IareaT_global = 0.0
19001903
do j=G%jsc,G%jec ; do i=G%isc,G%iec
1901-
tmpForSumming(i,j) = G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j)
1904+
tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j)
19021905
enddo ; enddo
19031906
G%areaT_global = reproducing_sum(tmpForSumming)
19041907
G%IareaT_global = 1. / (G%areaT_global)

‎src/parameterizations/lateral/MOM_MEKE.F90

+1-1
Original file line numberDiff line numberDiff line change
@@ -188,7 +188,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
188188
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)
189189
if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, scale=US%L_T_to_m_s**2)
190190
call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, scale=US%s_to_T)
191-
call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m)
191+
call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m*US%L_to_m**2)
192192
endif
193193

194194
sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping

‎src/parameterizations/vertical/MOM_set_diffusivity.F90

+2-2
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ module MOM_set_diffusivity
33

44
! This file is part of MOM6. See LICENSE.md for the license.
55

6+
use MOM_checksums, only : hchksum_pair
67
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
78
use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
89
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, &
344345

345346
if (CS%useKappaShear) then
346347
if (CS%debug) then
347-
call hchksum(u_h, "before calc_KS u_h",G%HI)
348-
call hchksum(v_h, "before calc_KS v_h",G%HI)
348+
call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, G%HI, scale=US%L_T_to_m_s)
349349
endif
350350
call cpu_clock_begin(id_clock_kappaShear)
351351
if (CS%Vertex_shear) then

‎src/parameterizations/vertical/MOM_sponge.F90

+1-1
Original file line numberDiff line numberDiff line change
@@ -420,7 +420,7 @@ subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml)
420420
eta_anom(i,j) = e_D(i,j,k) - CS%Ref_eta_im(j,k)
421421
if (CS%Ref_eta_im(j,K) < -G%bathyT(i,j)) eta_anom(i,j) = 0.0
422422
enddo ; enddo
423-
call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G)
423+
call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G, tmp_scale=US%Z_to_m)
424424
enddo
425425

426426
if (CS%fldno > 0) allocate(fld_mean_anom(G%isd:G%ied,nz,CS%fldno))

‎src/tracer/MOM_tracer_advect.F90

+9-9
Original file line numberDiff line numberDiff line change
@@ -485,8 +485,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
485485
else
486486
uhh(I) = uhr(I,j,k)
487487
endif
488-
!ts2(I) = 0.5*(1.0 + uhh(I)/(hprev(i+1,j,k)+h_neglect))
489-
CFL(I) = - uhh(I)/(hprev(i+1,j,k)+h_neglect) ! CFL is positive
488+
!ts2(I) = 0.5*(1.0 + uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)))
489+
CFL(I) = - uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)) ! CFL is positive
490490
else
491491
hup = hprev(i,j,k) - G%areaT(i,j)*min_h
492492
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, &
497497
else
498498
uhh(I) = uhr(I,j,k)
499499
endif
500-
!ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect))
501-
CFL(I) = uhh(I)/(hprev(i,j,k)+h_neglect) ! CFL is positive
500+
!ts2(I) = 0.5*(1.0 - uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)))
501+
CFL(I) = uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive
502502
endif
503503
enddo
504504

@@ -573,7 +573,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
573573
! Original implementation of PLM
574574
!flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I))
575575
endif
576-
!ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect))
576+
!ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j)))
577577
enddo ; enddo
578578
endif ! usePPM
579579

@@ -856,8 +856,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
856856
else
857857
vhh(i,J) = vhr(i,J,k)
858858
endif
859-
!ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k)+h_neglect))
860-
CFL(i) = - vhh(i,J) / (hprev(i,j+1,k)+h_neglect) ! CFL is positive
859+
!ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1))
860+
CFL(i) = - vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) ! CFL is positive
861861
else
862862
hup = hprev(i,j,k) - G%areaT(i,j)*min_h
863863
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, &
868868
else
869869
vhh(i,J) = vhr(i,J,k)
870870
endif
871-
!ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k)+h_neglect))
872-
CFL(i) = vhh(i,J) / (hprev(i,j,k)+h_neglect) ! CFL is positive
871+
!ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)))
872+
CFL(i) = vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive
873873
endif
874874
enddo
875875

0 commit comments

Comments
 (0)