From 8b5e10aa7040940bd3f8356b51f612eac0465580 Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Tue, 16 Nov 2021 21:25:15 -0500 Subject: [PATCH 1/8] (+) porous topography implementation (#3) * reads in porous topography parameters from CHANNEL_LIST_FILE *new module to compute curve fit for porous topography *porous constraints used to modify continuity_PPM, CoriolisAdv, and Rayleigh bottom channel drag --- src/core/MOM.F90 | 46 ++++- src/core/MOM_CoriolisAdv.F90 | 20 ++- src/core/MOM_continuity.F90 | 7 +- src/core/MOM_continuity_PPM.F90 | 145 ++++++++------- src/core/MOM_dynamics_split_RK2.F90 | 48 +++-- src/core/MOM_dynamics_unsplit.F90 | 43 ++++- src/core/MOM_dynamics_unsplit_RK2.F90 | 41 ++++- src/core/MOM_grid.F90 | 21 +++ src/core/MOM_porous_barriers.F90 | 166 ++++++++++++++++++ src/core/MOM_transcribe_grid.F90 | 16 ++ src/core/MOM_variables.F90 | 10 ++ src/framework/MOM_dyn_horgrid.F90 | 22 +++ .../MOM_shared_initialization.F90 | 55 +++++- .../vertical/MOM_set_viscosity.F90 | 15 +- 14 files changed, 545 insertions(+), 110 deletions(-) create mode 100644 src/core/MOM_porous_barriers.F90 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 5a1f4cf348..186c972b5e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -127,7 +127,7 @@ module MOM use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init use MOM_unit_scaling, only : unit_scaling_end, fix_restart_unit_scaling use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state -use MOM_variables, only : thermo_var_ptrs, vertvisc_type +use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_ptrs use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state use MOM_variables, only : rotate_surface_state use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd @@ -136,6 +136,8 @@ module MOM use MOM_wave_interface, only : wave_parameters_CS, waves_end use MOM_wave_interface, only : Update_Stokes_Drift +use MOM_porous_barriers, only : porous_widths + ! ODA modules use MOM_oda_driver_mod, only : ODA_CS, oda, init_oda, oda_end use MOM_oda_driver_mod, only : set_prior_tracer, set_analysis_time, apply_oda_tracer_increments @@ -399,6 +401,15 @@ module MOM type(ODA_CS), pointer :: odaCS => NULL() !< a pointer to the control structure for handling !! ensemble model state vectors and data assimilation !! increments and priors + type(porous_barrier_ptrs) :: pbv !< porous barrier fractional cell metrics + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) & + :: por_face_areaU !< fractional open area of U-faces [nondim] + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) & + :: por_face_areaV !< fractional open area of V-faces [nondim] + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) & + :: por_layer_widthU !< fractional open width of U-faces [nondim] + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) & + :: por_layer_widthV !< fractional open width of V-faces [nondim] type(particles), pointer :: particles => NULL() ! m or 1/eta_to_m] G => CS%G ; GV => CS%GV ; US => CS%US ; IDs => CS%IDs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -1047,13 +1060,16 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call diag_update_remap_grids(CS%diag) endif + !update porous barrier fractional cell metrics + call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv) + ! The bottom boundary layer properties need to be recalculated. if (bbl_time_int > 0.0) then call enable_averages(bbl_time_int, & Time_local + real_to_time(US%T_to_s*(bbl_time_int-dt)), CS%diag) ! Calculate the BBL properties and store them inside visc (u,h). call cpu_clock_begin(id_clock_BBL_visc) - call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, US, CS%set_visc_CSp) + call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM)") call disable_averaging(CS%diag) @@ -1076,7 +1092,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, & - CS%MEKE, CS%thickness_diffuse_CSp, waves=waves) + CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, waves=waves) if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_split (step_MOM)") elseif (CS%do_dynamics) then ! ------------------------------------ not SPLIT @@ -1090,11 +1106,11 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%use_RK2) then call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv) else call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, Waves=Waves) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, Waves=Waves) endif if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)") @@ -1296,6 +1312,9 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & integer :: halo_sz ! The size of a halo where data must be valid. integer :: i, j, k, is, ie, js, je, nz + real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_por ! layer interface heights + !! for porous topo. [Z ~> m or 1/eta_to_m] + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("step_MOM_thermo(), MOM.F90") @@ -1331,7 +1350,9 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & ! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics ! and set_viscous_BBL is called as a part of the dynamic stepping. call cpu_clock_begin(id_clock_BBL_visc) - call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp) + !update porous barrier fractional cell metrics + call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv) + call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM_thermo)") endif @@ -2330,6 +2351,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ALLOC_(CS%eta_av_bc(isd:ied,jsd:jed)) ; CS%eta_av_bc(:,:) = 0.0 ! -G%Z_ref CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0 + !allocate porous topography variables + ALLOC_(CS%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%por_face_areaU(:,:,:) = 1.0 + ALLOC_(CS%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%por_face_areaV(:,:,:) = 1.0 + ALLOC_(CS%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%por_layer_widthU(:,:,:) = 1.0 + ALLOC_(CS%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%por_layer_widthV(:,:,:) = 1.0 + CS%pbv%por_face_areaU => CS%por_face_areaU; CS%pbv%por_face_areaV=> CS%por_face_areaV + CS%pbv%por_layer_widthU => CS%por_layer_widthU; CS%pbv%por_layer_widthV => CS%por_layer_widthV ! Use the Wright equation of state by default, unless otherwise specified ! Note: this line and the following block ought to be in a separate ! initialization routine for tv. @@ -2647,7 +2675,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, & CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & - CS%visc, dirs, CS%ntrunc, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) + CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) if (CS%dtbt_reset_period > 0.0) then CS%dtbt_reset_interval = real_to_time(CS%dtbt_reset_period) ! Set dtbt_reset_time to be the next even multiple of dtbt_reset_interval. @@ -3581,6 +3609,10 @@ subroutine MOM_end(CS) if (CS%use_ALE_algorithm) call ALE_end(CS%ALE_CSp) + !deallocate porous topography variables + DEALLOC_(CS%por_face_areaU) ; DEALLOC_(CS%por_face_areaV) + DEALLOC_(CS%por_layer_widthU) ; DEALLOC_(CS%por_layer_widthV) + ! NOTE: Allocated in PressureForce_FV_Bouss if (associated(CS%tv%varT)) deallocate(CS%tv%varT) diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index e4d97ab53a..0de67d7159 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -14,7 +14,7 @@ module MOM_CoriolisAdv use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : accel_diag_ptrs +use MOM_variables, only : accel_diag_ptrs, porous_barrier_ptrs use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -117,7 +117,7 @@ module MOM_CoriolisAdv contains !> Calculates the Coriolis and momentum advection contributions to the acceleration. -subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) +subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, pbv) type(ocean_grid_type), intent(in) :: G !< Ocen grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1] @@ -135,6 +135,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) type(accel_diag_ptrs), intent(inout) :: AD !< Storage for acceleration diagnostics type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(CoriolisAdv_CS), pointer :: CS !< Control structure for MOM_CoriolisAdv + type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics ! Local variables real, dimension(SZIB_(G),SZJB_(G)) :: & @@ -285,7 +286,8 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) enddo ; enddo !$OMP parallel do default(private) shared(u,v,h,uh,vh,CAu,CAv,G,GV,CS,AD,Area_h,Area_q,& - !$OMP RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,vol_neglect,h_tiny,OBC,eps_vel) + !$OMP RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,vol_neglect,h_tiny,OBC,eps_vel, & + !$OMP pbv) do k=1,nz ! Here the second order accurate layer potential vorticities, q, @@ -306,10 +308,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) if (CS%Coriolis_En_Dis) then do j=Jsq,Jeq+1 ; do I=is-1,ie - uh_center(I,j) = 0.5 * (G%dy_Cu(I,j) * u(I,j,k)) * (h(i,j,k) + h(i+1,j,k)) + uh_center(I,j) = 0.5 * ((G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k)) * u(I,j,k)) * (h(i,j,k) + h(i+1,j,k)) enddo ; enddo do J=js-1,je ; do i=Isq,Ieq+1 - vh_center(i,J) = 0.5 * (G%dx_Cv(i,J) * v(i,J,k)) * (h(i,j,k) + h(i,j+1,k)) + vh_center(i,J) = 0.5 * ((G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k)) * v(i,J,k)) * (h(i,j,k) + h(i,j+1,k)) enddo ; enddo endif @@ -352,9 +354,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) if (CS%Coriolis_En_Dis) then do i = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+2,OBC%segment(n)%HI%ied) if (OBC%segment(n)%direction == OBC_DIRECTION_N) then - vh_center(i,J) = G%dx_Cv(i,J) * v(i,J,k) * h(i,j,k) + vh_center(i,J) = (G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k)) * v(i,J,k) * h(i,j,k) else ! (OBC%segment(n)%direction == OBC_DIRECTION_S) - vh_center(i,J) = G%dx_Cv(i,J) * v(i,J,k) * h(i,j+1,k) + vh_center(i,J) = (G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k)) * v(i,J,k) * h(i,j+1,k) endif enddo endif @@ -391,9 +393,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) if (CS%Coriolis_En_Dis) then do j = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+2,OBC%segment(n)%HI%jed) if (OBC%segment(n)%direction == OBC_DIRECTION_E) then - uh_center(I,j) = G%dy_Cu(I,j) * u(I,j,k) * h(i,j,k) + uh_center(I,j) = (G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k)) * u(I,j,k) * h(i,j,k) else ! (OBC%segment(n)%direction == OBC_DIRECTION_W) - uh_center(I,j) = G%dy_Cu(I,j) * u(I,j,k) * h(i+1,j,k) + uh_center(I,j) = (G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k)) * u(I,j,k) * h(i+1,j,k) endif enddo endif diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 655055b03d..20f9c34d65 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -13,7 +13,7 @@ module MOM_continuity use MOM_grid, only : ocean_grid_type use MOM_open_boundary, only : ocean_OBC_type use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : BT_cont_type +use MOM_variables, only : BT_cont_type, porous_barrier_ptrs use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -39,7 +39,7 @@ module MOM_continuity !> Time steps the layer thicknesses, using a monotonically limited, directionally split PPM scheme, !! based on Lin (1994). -subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt, & +subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, & visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. @@ -61,6 +61,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt, type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_CS), pointer :: CS !< Control structure for mom_continuity. type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics real, dimension(SZIB_(G),SZJ_(G)), & optional, intent(in) :: uhbt !< The vertically summed volume !! flux through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -95,7 +96,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt, " one must be present in call to continuity.") if (CS%continuity_scheme == PPM_SCHEME) then - call continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS%PPM_CSp, OBC, uhbt, vhbt, & + call continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS%PPM_CSp, OBC, pbv, uhbt, vhbt, & visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont=BT_cont) else call MOM_error(FATAL, "continuity: Unrecognized value of continuity_scheme") diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index d30e1af0f2..c4c96fa392 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -11,7 +11,7 @@ module MOM_continuity_PPM use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : BT_cont_type +use MOM_variables, only : BT_cont_type, porous_barrier_ptrs use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -73,7 +73,7 @@ module MOM_continuity_PPM !> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme, !! based on Lin (1994). -subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt, & +subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, & visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. @@ -93,6 +93,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_PPM_CS), pointer :: CS !< Module's control structure. type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + type(porous_barrier_ptrs), intent(in) :: pbv !< pointers to porous barrier fractional cell metrics real, dimension(SZIB_(G),SZJ_(G)), & optional, intent(in) :: uhbt !< The summed volume flux through zonal faces !! [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -148,7 +149,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh ! First, advect zonally. LB%ish = G%isc ; LB%ieh = G%iec LB%jsh = G%jsc-stencil ; LB%jeh = G%jec+stencil - call zonal_mass_flux(u, hin, uh, dt, G, GV, US, CS, LB, OBC, uhbt, visc_rem_u, u_cor, BT_cont) + call zonal_mass_flux(u, hin, uh, dt, G, GV, US, CS, LB, OBC, & + pbv%por_face_areaU, uhbt, visc_rem_u, u_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -163,7 +165,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh ! Now advect meridionally, using the updated thicknesses to determine ! the fluxes. - call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS, LB, OBC, vhbt, visc_rem_v, v_cor, BT_cont) + call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS, LB, OBC, & + pbv%por_face_areaV, vhbt, visc_rem_v, v_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -179,7 +182,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh LB%ish = G%isc-stencil ; LB%ieh = G%iec+stencil LB%jsh = G%jsc ; LB%jeh = G%jec - call meridional_mass_flux(v, hin, vh, dt, G, GV, US, CS, LB, OBC, vhbt, visc_rem_v, v_cor, BT_cont) + call meridional_mass_flux(v, hin, vh, dt, G, GV, US, CS, LB, OBC, & + pbv%por_face_areaV, vhbt, visc_rem_v, v_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -191,7 +195,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh ! Now advect zonally, using the updated thicknesses to determine ! the fluxes. LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec - call zonal_mass_flux(u, h, uh, dt, G, GV, US, CS, LB, OBC, uhbt, visc_rem_u, u_cor, BT_cont) + call zonal_mass_flux(u, h, uh, dt, G, GV, US, CS, LB, OBC, & + pbv%por_face_areaU, uhbt, visc_rem_u, u_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -207,7 +212,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh end subroutine continuity_PPM !> Calculates the mass or volume fluxes through the zonal faces, and other related quantities. -subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & +subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_areaU, uhbt, & visc_rem_u, u_cor, BT_cont) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. @@ -223,6 +228,8 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & type(continuity_PPM_CS), pointer :: CS !< This module's control structure. type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), & + intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim] real, dimension(SZIB_(G),SZJ_(G)), & optional, intent(in) :: uhbt !< The summed volume flux through zonal faces !! [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -300,7 +307,8 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & call cpu_clock_begin(id_clock_correct) !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,u,h_in,h_L,h_R,use_visc_rem,visc_rem_u, & !$OMP uh,dt,US,G,GV,CS,local_specified_BC,OBC,uhbt,set_BT_cont, & -!$OMP CFL_dt,I_dt,u_cor,BT_cont, local_Flather_OBC) & +!$OMP CFL_dt,I_dt,u_cor,BT_cont, local_Flather_OBC, & +!$OMP por_face_areaU) & !$OMP private(do_I,duhdu,du,du_max_CFL,du_min_CFL,uh_tot_0,duhdu_tot_0, & !$OMP is_simple,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W, & !$OMP any_simple_OBC,l_seg) & @@ -315,7 +323,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & enddo ; endif call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_L(:,j,k), h_R(:,j,k), & uh(:,j,k), duhdu(:,k), visc_rem(:,k), & - dt, G, US, j, ish, ieh, do_I, CS%vol_CFL, OBC) + dt, G, US, j, ish, ieh, do_I, CS%vol_CFL, por_face_areaU(:,j,k), OBC) if (local_specified_BC) then do I=ish-1,ieh l_seg = OBC%segnum_u(I,j) @@ -428,7 +436,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & if (present(uhbt)) then call zonal_flux_adjust(u, h_in, h_L, h_R, uhbt(:,j), uh_tot_0, duhdu_tot_0, du, & du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I, uh, OBC=OBC) + j, ish, ieh, do_I, por_face_areaU, uh, OBC=OBC) if (present(u_cor)) then ; do k=1,nz do I=ish-1,ieh ; u_cor(I,j,k) = u(I,j,k) + du(I) * visc_rem(I,k) ; enddo @@ -447,7 +455,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & if (set_BT_cont) then call set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0,& du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - visc_rem_max, j, ish, ieh, do_I) + visc_rem_max, j, ish, ieh, do_I, por_face_areaU) if (any_simple_OBC) then do I=ish-1,ieh l_seg = OBC%segnum_u(I,j) @@ -456,7 +464,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & if (l_seg /= OBC_NONE) & do_I(I) = OBC%segment(l_seg)%specified - if (do_I(I)) FAuI(I) = GV%H_subroundoff*G%dy_Cu(I,j) + if (do_I(I)) FAuI(I) = GV%H_subroundoff*(G%dy_Cu(I,j)*por_face_areaU(I,j,k)) enddo ! NOTE: do_I(I) should prevent access to segment OBC_NONE do k=1,nz ; do I=ish-1,ieh ; if (do_I(I)) then @@ -484,7 +492,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & if (OBC%segment(n)%direction == OBC_DIRECTION_E) then do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed FA_u = 0.0 - do k=1,nz ; FA_u = FA_u + h_in(i,j,k)*G%dy_Cu(I,j) ; enddo + do k=1,nz ; FA_u = FA_u + h_in(i,j,k)*(G%dy_Cu(I,j)*por_face_areaU(I,j,k)) ; enddo BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u BT_cont%FA_u_WW(I,j) = FA_u ; BT_cont%FA_u_EE(I,j) = FA_u BT_cont%uBT_WW(I,j) = 0.0 ; BT_cont%uBT_EE(I,j) = 0.0 @@ -492,7 +500,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & else do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed FA_u = 0.0 - do k=1,nz ; FA_u = FA_u + h_in(i+1,j,k)*G%dy_Cu(I,j) ; enddo + do k=1,nz ; FA_u = FA_u + h_in(i+1,j,k)*(G%dy_Cu(I,j)*por_face_areaU(I,j,k)) ; enddo BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u BT_cont%FA_u_WW(I,j) = FA_u ; BT_cont%FA_u_EE(I,j) = FA_u BT_cont%uBT_WW(I,j) = 0.0 ; BT_cont%uBT_EE(I,j) = 0.0 @@ -506,10 +514,10 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & if (set_BT_cont) then ; if (allocated(BT_cont%h_u)) then if (present(u_cor)) then call zonal_face_thickness(u_cor, h_in, h_L, h_R, BT_cont%h_u, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_u) + CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU(:,j,k), visc_rem_u) else call zonal_face_thickness(u, h_in, h_L, h_R, BT_cont%h_u, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_u) + CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaU(:,j,k), visc_rem_u) endif endif ; endif @@ -517,7 +525,7 @@ end subroutine zonal_mass_flux !> Evaluates the zonal mass or volume fluxes in a layer. subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, & - ish, ieh, do_I, vol_CFL, OBC) + ish, ieh, do_I, vol_CFL, por_face_areaU, OBC) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. real, dimension(SZIB_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. real, dimension(SZIB_(G)), intent(in) :: visc_rem !< Both the fraction of the @@ -539,6 +547,7 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, & integer, intent(in) :: ieh !< End of index range. logical, dimension(SZIB_(G)), intent(in) :: do_I !< Which i values to work on. logical, intent(in) :: vol_CFL !< If true, rescale the + real, dimension(SZIB_(G)), intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim] !! ratio of face areas to the cell areas when estimating the CFL number. type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. ! Local variables @@ -561,21 +570,21 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, & if (vol_CFL) then ; CFL = (u(I) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j)) else ; CFL = u(I) * dt * G%IdxT(i,j) ; endif curv_3 = h_L(i) + h_R(i) - 2.0*h(i) - uh(I) = G%dy_Cu(I,j) * u(I) * & + uh(I) = (G%dy_Cu(I,j) * por_face_areaU(I))* u(I) * & (h_R(i) + CFL * (0.5*(h_L(i) - h_R(i)) + curv_3*(CFL - 1.5))) h_marg = h_R(i) + CFL * ((h_L(i) - h_R(i)) + 3.0*curv_3*(CFL - 1.0)) elseif (u(I) < 0.0) then if (vol_CFL) then ; CFL = (-u(I) * dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j)) else ; CFL = -u(I) * dt * G%IdxT(i+1,j) ; endif curv_3 = h_L(i+1) + h_R(i+1) - 2.0*h(i+1) - uh(I) = G%dy_Cu(I,j) * u(I) * & + uh(I) = (G%dy_Cu(I,j) * por_face_areaU(I)) * u(I) * & (h_L(i+1) + CFL * (0.5*(h_R(i+1)-h_L(i+1)) + curv_3*(CFL - 1.5))) h_marg = h_L(i+1) + CFL * ((h_R(i+1)-h_L(i+1)) + 3.0*curv_3*(CFL - 1.0)) else uh(I) = 0.0 h_marg = 0.5 * (h_L(i+1) + h_R(i)) endif - duhdu(I) = G%dy_Cu(I,j) * h_marg * visc_rem(I) + duhdu(I) = (G%dy_Cu(I,j)* por_face_areaU(I)) * h_marg * visc_rem(I) endif ; enddo if (local_open_BC) then @@ -585,11 +594,11 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, & if (l_seg /= OBC_NONE) then if (OBC%segment(l_seg)%open) then if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then - uh(I) = G%dy_Cu(I,j) * u(I) * h(i) - duhdu(I) = G%dy_Cu(I,j) * h(i) * visc_rem(I) + uh(I) = (G%dy_Cu(I,j)* por_face_areaU(I)) * u(I) * h(i) + duhdu(I) = (G%dy_Cu(I,j)* por_face_areaU(I)) * h(i) * visc_rem(I) else - uh(I) = G%dy_Cu(I,j) * u(I) * h(i+1) - duhdu(I) = G%dy_Cu(I,j) * h(i+1) * visc_rem(I) + uh(I) = (G%dy_Cu(I,j)* por_face_areaU(I)) * u(I) * h(i+1) + duhdu(I) = (G%dy_Cu(I,j)* por_face_areaU(I)) * h(i+1) * visc_rem(I) endif endif endif @@ -599,7 +608,7 @@ end subroutine zonal_flux_layer !> Sets the effective interface thickness at each zonal velocity point. subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, & - marginal, OBC, visc_rem_u) + marginal, OBC, por_face_areaU, visc_rem_u) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. @@ -617,6 +626,8 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, !! of face areas to the cell areas when estimating the CFL number. logical, intent(in) :: marginal !< If true, report the !! marginal face thicknesses; otherwise report transport-averaged thicknesses. + real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), & + intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim] type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: visc_rem_u @@ -706,7 +717,8 @@ end subroutine zonal_face_thickness !! desired barotropic (layer-summed) transport. subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & du, du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I_in, uh_3d, OBC) + j, ish, ieh, do_I_in, por_face_areaU, uh_3d, OBC) + type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. @@ -742,6 +754,8 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & integer, intent(in) :: ieh !< End of index range. logical, dimension(SZIB_(G)), intent(in) :: do_I_in !< !! A logical flag indicating which I values to work on. + real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), & + intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), optional, intent(inout) :: uh_3d !< !! Volume flux through zonal faces = u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]. type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -826,7 +840,7 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & do I=ish-1,ieh ; u_new(I) = u(I,j,k) + du(I) * visc_rem(I,k) ; enddo call zonal_flux_layer(u_new, h_in(:,j,k), h_L(:,j,k), h_R(:,j,k), & uh_aux(:,k), duhdu(:,k), visc_rem(:,k), & - dt, G, US, j, ish, ieh, do_I, CS%vol_CFL, OBC) + dt, G, US, j, ish, ieh, do_I, CS%vol_CFL,por_face_areaU(:,j,k),OBC) enddo ; endif if (itt < max_itts) then @@ -856,7 +870,7 @@ end subroutine zonal_flux_adjust !! function of barotropic flow to agree closely with the sum of the layer's transports. subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, & du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - visc_rem_max, j, ish, ieh, do_I) + visc_rem_max, j, ish, ieh, do_I, por_face_areaU) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. @@ -890,6 +904,8 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, integer, intent(in) :: ieh !< End of index range. logical, dimension(SZIB_(G)), intent(in) :: do_I !< A logical flag indicating !! which I values to work on. + real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), & + intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim] ! Local variables real, dimension(SZIB_(G)) :: & du0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1]. @@ -931,7 +947,7 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, do I=ish-1,ieh ; zeros(I) = 0.0 ; enddo call zonal_flux_adjust(u, h_in, h_L, h_R, zeros, uh_tot_0, duhdu_tot_0, du0, & du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I) + j, ish, ieh, do_I, por_face_areaU) ! Determine the westerly- and easterly- fluxes. Choose a sufficiently ! negative velocity correction for the easterly-flux, and a sufficiently @@ -972,11 +988,11 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, u_0(I) = u(I,j,k) + du0(I) * visc_rem(I,k) endif ; enddo call zonal_flux_layer(u_0, h_in(:,j,k), h_L(:,j,k), h_R(:,j,k), uh_0, duhdu_0, & - visc_rem(:,k), dt, G, US, j, ish, ieh, do_I, CS%vol_CFL) + visc_rem(:,k), dt, G, US, j, ish, ieh, do_I, CS%vol_CFL, por_face_areaU(:,j,k)) call zonal_flux_layer(u_L, h_in(:,j,k), h_L(:,j,k), h_R(:,j,k), uh_L, duhdu_L, & - visc_rem(:,k), dt, G, US, j, ish, ieh, do_I, CS%vol_CFL) + visc_rem(:,k), dt, G, US, j, ish, ieh, do_I, CS%vol_CFL,por_face_areaU(:,j,k)) call zonal_flux_layer(u_R, h_in(:,j,k), h_L(:,j,k), h_R(:,j,k), uh_R, duhdu_R, & - visc_rem(:,k), dt, G, US, j, ish, ieh, do_I, CS%vol_CFL) + visc_rem(:,k), dt, G, US, j, ish, ieh, do_I, CS%vol_CFL,por_face_areaU(:,j,k)) do I=ish-1,ieh ; if (do_I(I)) then FAmt_0(I) = FAmt_0(I) + duhdu_0(I) FAmt_L(I) = FAmt_L(I) + duhdu_L(I) @@ -1018,7 +1034,7 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, end subroutine set_zonal_BT_cont !> Calculates the mass or volume fluxes through the meridional faces, and other related quantities. -subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & +subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_face_areaV, vhbt, & visc_rem_v, v_cor, BT_cont) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. @@ -1026,13 +1042,15 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to !! calculate fluxes [H ~> m or kg m-2] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: vh !< Volume flux through meridional - !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1] + !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1]. real, intent(in) :: dt !< Time increment [T ~> s]. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_PPM_CS), pointer :: CS !< This module's control structure.G type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. type(ocean_OBC_type), pointer :: OBC !< Open boundary condition type !! specifies whether, where, and what open boundary conditions are used. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim] real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vhbt !< The summed volume flux through !< meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -1110,7 +1128,8 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & call cpu_clock_begin(id_clock_correct) !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,v,h_in,h_L,h_R,vh,use_visc_rem, & !$OMP visc_rem_v,dt,US,G,GV,CS,local_specified_BC,OBC,vhbt, & -!$OMP set_BT_cont,CFL_dt,I_dt,v_cor,BT_cont, local_Flather_OBC ) & +!$OMP set_BT_cont,CFL_dt,I_dt,v_cor,BT_cont, local_Flather_OBC, & +!$OMP por_face_areaV) & !$OMP private(do_I,dvhdv,dv,dv_max_CFL,dv_min_CFL,vh_tot_0, & !$OMP dvhdv_tot_0,visc_rem_max,I_vrm,dv_lim,dy_N, & !$OMP is_simple,FAvi,dy_S,any_simple_OBC,l_seg) & @@ -1125,7 +1144,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & enddo ; endif call merid_flux_layer(v(:,J,k), h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), & vh(:,J,k), dvhdv(:,k), visc_rem(:,k), & - dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, OBC) + dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, por_face_areaV(:,:,k), OBC) if (local_specified_BC) then do i=ish,ieh l_seg = OBC%segnum_v(i,J) @@ -1234,7 +1253,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & if (present(vhbt)) then call meridional_flux_adjust(v, h_in, h_L, h_R, vhbt(:,J), vh_tot_0, dvhdv_tot_0, dv, & dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I, vh, OBC=OBC) + j, ish, ieh, do_I, por_face_areaV, vh, OBC=OBC) if (present(v_cor)) then ; do k=1,nz do i=ish,ieh ; v_cor(i,J,k) = v(i,J,k) + dv(i) * visc_rem(i,k) ; enddo @@ -1252,7 +1271,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & if (set_BT_cont) then call set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0,& dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - visc_rem_max, J, ish, ieh, do_I) + visc_rem_max, J, ish, ieh, do_I, por_face_areaV) if (any_simple_OBC) then do i=ish,ieh l_seg = OBC%segnum_v(i,J) @@ -1261,7 +1280,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & if(l_seg /= OBC_NONE) & do_I(i) = (OBC%segment(l_seg)%specified) - if (do_I(i)) FAvi(i) = GV%H_subroundoff*G%dx_Cv(i,J) + if (do_I(i)) FAvi(i) = GV%H_subroundoff*(G%dx_Cv(i,J)*por_face_areaV(i,J,k)) enddo ! NOTE: do_I(I) should prevent access to segment OBC_NONE do k=1,nz ; do i=ish,ieh ; if (do_I(i)) then @@ -1289,7 +1308,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & if (OBC%segment(n)%direction == OBC_DIRECTION_N) then do i = OBC%segment(n)%HI%Isd, OBC%segment(n)%HI%Ied FA_v = 0.0 - do k=1,nz ; FA_v = FA_v + h_in(i,j,k)*G%dx_Cv(i,J) ; enddo + do k=1,nz ; FA_v = FA_v + h_in(i,j,k)*(G%dx_Cv(i,J)*por_face_areaV(i,J,k)) ; enddo BT_cont%FA_v_S0(i,J) = FA_v ; BT_cont%FA_v_N0(i,J) = FA_v BT_cont%FA_v_SS(i,J) = FA_v ; BT_cont%FA_v_NN(i,J) = FA_v BT_cont%vBT_SS(i,J) = 0.0 ; BT_cont%vBT_NN(i,J) = 0.0 @@ -1297,7 +1316,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & else do i = OBC%segment(n)%HI%Isd, OBC%segment(n)%HI%Ied FA_v = 0.0 - do k=1,nz ; FA_v = FA_v + h_in(i,j+1,k)*G%dx_Cv(i,J) ; enddo + do k=1,nz ; FA_v = FA_v + h_in(i,j+1,k)*(G%dx_Cv(i,J)*por_face_areaV(i,J,k)) ; enddo BT_cont%FA_v_S0(i,J) = FA_v ; BT_cont%FA_v_N0(i,J) = FA_v BT_cont%FA_v_SS(i,J) = FA_v ; BT_cont%FA_v_NN(i,J) = FA_v BT_cont%vBT_SS(i,J) = 0.0 ; BT_cont%vBT_NN(i,J) = 0.0 @@ -1311,10 +1330,10 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & if (set_BT_cont) then ; if (allocated(BT_cont%h_v)) then if (present(v_cor)) then call merid_face_thickness(v_cor, h_in, h_L, h_R, BT_cont%h_v, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_v) + CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaV, visc_rem_v) else call merid_face_thickness(v, h_in, h_L, h_R, BT_cont%h_v, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_v) + CS%vol_CFL, CS%marginal_faces, OBC, por_face_areaV, visc_rem_v) endif endif ; endif @@ -1322,7 +1341,7 @@ end subroutine meridional_mass_flux !> Evaluates the meridional mass or volume fluxes in a layer. subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & - ish, ieh, do_I, vol_CFL, OBC) + ish, ieh, do_I, vol_CFL, por_face_areaV, OBC) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. real, dimension(SZI_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G)), intent(in) :: visc_rem !< Both the fraction of the @@ -1348,6 +1367,8 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & logical, dimension(SZI_(G)), intent(in) :: do_I !< Which i values to work on. logical, intent(in) :: vol_CFL !< If true, rescale the !! ratio of face areas to the cell areas when estimating the CFL number. + real, dimension(SZI_(G), SZJB_(G)), & + intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim] type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim] @@ -1368,7 +1389,7 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & if (vol_CFL) then ; CFL = (v(i) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j)) else ; CFL = v(i) * dt * G%IdyT(i,j) ; endif curv_3 = h_L(i,j) + h_R(i,j) - 2.0*h(i,j) - vh(i) = G%dx_Cv(i,J) * v(i) * ( h_R(i,j) + CFL * & + vh(i) = (G%dx_Cv(i,J)*por_face_areaV(i,J)) * v(i) * ( h_R(i,j) + CFL * & (0.5*(h_L(i,j) - h_R(i,j)) + curv_3*(CFL - 1.5)) ) h_marg = h_R(i,j) + CFL * ((h_L(i,j) - h_R(i,j)) + & 3.0*curv_3*(CFL - 1.0)) @@ -1376,7 +1397,7 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & if (vol_CFL) then ; CFL = (-v(i) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1)) else ; CFL = -v(i) * dt * G%IdyT(i,j+1) ; endif curv_3 = h_L(i,j+1) + h_R(i,j+1) - 2.0*h(i,j+1) - vh(i) = G%dx_Cv(i,J) * v(i) * ( h_L(i,j+1) + CFL * & + vh(i) = (G%dx_Cv(i,J)*por_face_areaV(i,J)) * v(i) * ( h_L(i,j+1) + CFL * & (0.5*(h_R(i,j+1)-h_L(i,j+1)) + curv_3*(CFL - 1.5)) ) h_marg = h_L(i,j+1) + CFL * ((h_R(i,j+1)-h_L(i,j+1)) + & 3.0*curv_3*(CFL - 1.0)) @@ -1384,7 +1405,7 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & vh(i) = 0.0 h_marg = 0.5 * (h_L(i,j+1) + h_R(i,j)) endif - dvhdv(i) = G%dx_Cv(i,J) * h_marg * visc_rem(i) + dvhdv(i) = (G%dx_Cv(i,J)*por_face_areaV(i,J)) * h_marg * visc_rem(i) endif ; enddo if (local_open_BC) then @@ -1394,11 +1415,11 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & if (l_seg /= OBC_NONE) then if (OBC%segment(l_seg)%open) then if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then - vh(i) = G%dx_Cv(i,J) * v(i) * h(i,j) - dvhdv(i) = G%dx_Cv(i,J) * h(i,j) * visc_rem(i) + vh(i) = (G%dx_Cv(i,J)*por_face_areaV(i,J)) * v(i) * h(i,j) + dvhdv(i) = (G%dx_Cv(i,J)*por_face_areaV(i,J)) * h(i,j) * visc_rem(i) else - vh(i) = G%dx_Cv(i,J) * v(i) * h(i,j+1) - dvhdv(i) = G%dx_Cv(i,J) * h(i,j+1) * visc_rem(i) + vh(i) = (G%dx_Cv(i,J)*por_face_areaV(i,J)) * v(i) * h(i,j+1) + dvhdv(i) = (G%dx_Cv(i,J)*por_face_areaV(i,J)) * h(i,j+1) * visc_rem(i) endif endif endif @@ -1408,7 +1429,7 @@ end subroutine merid_flux_layer !> Sets the effective interface thickness at each meridional velocity point. subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, & - marginal, OBC, visc_rem_v) + marginal, OBC, por_face_areaV, visc_rem_v) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]. @@ -1428,6 +1449,8 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, logical, intent(in) :: marginal !< If true, report the marginal !! face thicknesses; otherwise report transport-averaged thicknesses. type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), optional, intent(in) :: visc_rem_v !< Both the fraction !! of the momentum originally in a layer that remains after a time-step of !! viscosity, and the fraction of a time-step's worth of a barotropic @@ -1516,7 +1539,7 @@ end subroutine merid_face_thickness !> Returns the barotropic velocity adjustment that gives the desired barotropic (layer-summed) transport. subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0, & dv, dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I_in, vh_3d, OBC) + j, ish, ieh, do_I_in, por_face_areaV, vh_3d, OBC) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -1551,6 +1574,8 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 integer, intent(in) :: ieh !< End of index range. logical, dimension(SZI_(G)), & intent(in) :: do_I_in !< A flag indicating which I values to work on. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & optional, intent(inout) :: vh_3d !< Volume flux through meridional !! faces = v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -1636,7 +1661,7 @@ subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0 do i=ish,ieh ; v_new(i) = v(i,J,k) + dv(i) * visc_rem(i,k) ; enddo call merid_flux_layer(v_new, h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), & vh_aux(:,k), dvhdv(:,k), visc_rem(:,k), & - dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, OBC) + dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, por_face_areaV(:,:,k), OBC) enddo ; endif if (itt < max_itts) then @@ -1666,7 +1691,7 @@ end subroutine meridional_flux_adjust !! function of barotropic flow to agree closely with the sum of the layer's transports. subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, & dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - visc_rem_max, j, ish, ieh, do_I) + visc_rem_max, j, ish, ieh, do_I, por_face_areaV) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]. @@ -1700,6 +1725,8 @@ subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, integer, intent(in) :: ieh !< End of index range. logical, dimension(SZI_(G)), intent(in) :: do_I !< A logical flag indicating !! which I values to work on. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim] ! Local variables real, dimension(SZI_(G)) :: & dv0, & ! The barotropic velocity increment that gives 0 transport [L T-1 ~> m s-1]. @@ -1741,7 +1768,7 @@ subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, do i=ish,ieh ; zeros(i) = 0.0 ; enddo call meridional_flux_adjust(v, h_in, h_L, h_R, zeros, vh_tot_0, dvhdv_tot_0, dv0, & dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & - j, ish, ieh, do_I) + j, ish, ieh, do_I, por_face_areaV) ! Determine the southerly- and northerly- fluxes. Choose a sufficiently ! negative velocity correction for the northerly-flux, and a sufficiently @@ -1782,11 +1809,11 @@ subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, v_0(i) = v(I,j,k) + dv0(i) * visc_rem(i,k) endif ; enddo call merid_flux_layer(v_0, h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), vh_0, dvhdv_0, & - visc_rem(:,k), dt, G, US, J, ish, ieh, do_I, CS%vol_CFL) + visc_rem(:,k), dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, por_face_areaV(:,:,k)) call merid_flux_layer(v_L, h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), vh_L, dvhdv_L, & - visc_rem(:,k), dt, G, US, J, ish, ieh, do_I, CS%vol_CFL) + visc_rem(:,k), dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, por_face_areaV(:,:,k)) call merid_flux_layer(v_R, h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), vh_R, dvhdv_R, & - visc_rem(:,k), dt, G, US, J, ish, ieh, do_I, CS%vol_CFL) + visc_rem(:,k), dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, por_face_areaV(:,:,k)) do i=ish,ieh ; if (do_I(i)) then FAmt_0(i) = FAmt_0(i) + dvhdv_0(i) FAmt_L(i) = FAmt_L(i) + dvhdv_L(i) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index f9d70d65d7..fa796c7fe8 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -3,7 +3,7 @@ module MOM_dynamics_split_RK2 ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_variables, only : vertvisc_type, thermo_var_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_ptrs use MOM_variables, only : BT_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs use MOM_forcing_type, only : mech_forcing @@ -165,6 +165,7 @@ module MOM_dynamics_split_RK2 integer :: id_umo_2d = -1, id_vmo_2d = -1 integer :: id_PFu = -1, id_PFv = -1 integer :: id_CAu = -1, id_CAv = -1 + integer :: id_ueffA = -1, id_veffA = -1 ! integer :: id_hf_PFu = -1, id_hf_PFv = -1 integer :: id_h_PFu = -1, id_h_PFv = -1 integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1 @@ -255,7 +256,7 @@ module MOM_dynamics_split_RK2 !> RK2 splitting for time stepping MOM adiabatic dynamics subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, & uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, calc_dtbt, VarMix, & - MEKE, thickness_diffuse_CSp, Waves) + MEKE, thickness_diffuse_CSp, pbv, Waves) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -294,6 +295,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s type(MEKE_type), pointer :: MEKE !< related to mesoscale eddy kinetic energy param type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp !< Pointer to a structure containing !! interface height diffusivities + type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -301,7 +303,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted thickness [H ~> m or kg m-2]. - + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffA ! Effective Area of U-Faces [H L ~> m2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffA ! Effective Area of V-Faces [H L ~> m2] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_bc_accel real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_bc_accel ! u_bc_accel and v_bc_accel are the summed baroclinic accelerations of each @@ -395,6 +398,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s do j=G%jsdB,G%jedB ; do i=G%isd,G%ied ; vp(i,j,k) = 0.0 ; enddo ; enddo do j=G%jsd,G%jed ; do i=G%isd,G%ied ; hp(i,j,k) = h(i,j,k) ; enddo ; enddo enddo + ueffA(:,:,:) = 0; veffA(:,:,:) = 0 ! Update CFL truncation value as function of time call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp) @@ -485,7 +489,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av call cpu_clock_begin(id_clock_Cor) call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, Gv, US, CS%CoriolisAdv_CSp) + G, Gv, US, CS%CoriolisAdv_CSp, pbv) call cpu_clock_end(id_clock_Cor) if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)") @@ -562,7 +566,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! u_accel_bt = layer accelerations due to barotropic solver if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, & + call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, & visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (BT_cont_BT_thick) then @@ -647,7 +651,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! uh = u_av * h ! hp = h + dt * div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, & + call continuity(up, vp, h, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, & CS%uhbt, CS%vhbt, CS%visc_rem_u, CS%visc_rem_v, & u_av, v_av, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) @@ -738,7 +742,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av call cpu_clock_begin(id_clock_Cor) call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, GV, US, CS%CoriolisAdv_CSp) + G, GV, US, CS%CoriolisAdv_CSp, pbv) call cpu_clock_end(id_clock_Cor) if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)") @@ -850,7 +854,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! h = h + dt * div . uh ! u_av and v_av adjusted so their mass transports match uhbt and vhbt. call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, & + call continuity(u, v, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, & CS%uhbt, CS%vhbt, CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) call do_group_pass(CS%pass_h, G%Domain, clock=id_clock_pass) @@ -906,6 +910,22 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag) if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag) + ! Calculate effective areas and post data + if (CS%id_ueffA > 0) then + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_ueffA, ueffA, CS%diag) + endif + + if (CS%id_veffA > 0) then + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_veffA, veffA, CS%diag) + endif + + ! Diagnostics for terms multiplied by fractional thicknesses ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option. @@ -1238,7 +1258,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, & - visc, dirs, ntrunc, calc_dtbt, cont_stencil) + visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil) type(ocean_grid_type), intent(inout) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1278,6 +1298,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !! the number of times the velocity is !! truncated (this should be 0). logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step + type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics integer, intent(out) :: cont_stencil !< The stencil for thickness !! from the continuity solver. @@ -1474,7 +1495,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param if (.not. query_initialized(uh,"uh",restart_CS) .or. & .not. query_initialized(vh,"vh",restart_CS)) then do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo - call continuity(u, v, h, h_tmp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC) + call continuity(u, v, h, h_tmp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv) call pass_var(h_tmp, G%Domain, clock=id_clock_pass_init) do k=1,nz ; do j=jsd,jed ; do i=isd,ied CS%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k)) @@ -1519,6 +1540,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 'Zonal Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & + 'Effective U-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + y_cell_method='sum', v_extensive = .true.) + CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & + 'Effective V-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + x_cell_method='sum', v_extensive = .true.) + !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, & ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', & diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 48d767e1a8..7cfc9d649c 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -50,7 +50,7 @@ module MOM_dynamics_unsplit !* * !********+*********+*********+*********+*********+*********+*********+** -use MOM_variables, only : vertvisc_type, thermo_var_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_ptrs use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs use MOM_forcing_type, only : mech_forcing use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum @@ -128,6 +128,7 @@ module MOM_dynamics_unsplit !>@{ Diagnostic IDs integer :: id_uh = -1, id_vh = -1 + integer :: id_ueffA = -1, id_veffA = -1 integer :: id_PFu = -1, id_PFv = -1, id_CAu = -1, id_CAv = -1 !>@} @@ -186,7 +187,7 @@ module MOM_dynamics_unsplit !! 3rd order (for the inviscid momentum equations) order scheme subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE, Waves) + VarMix, MEKE, pbv, Waves) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -221,6 +222,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & !! that specify the spatially variable viscosities. type(MEKE_type), pointer :: MEKE !< A pointer to a structure containing !! fields related to the Mesoscale Eddy Kinetic Energy. + type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -228,6 +230,8 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av, hp ! Prediced or averaged layer thicknesses [H ~> m or kg m-2] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up, upp ! Predicted zonal velocities [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp, vpp ! Predicted meridional velocities [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffA ! Effective Area of U-Faces [H L ~> m2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffA ! Effective Area of V-Faces [H L ~> m2] real, dimension(:,:), pointer :: p_surf => NULL() real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. real :: dt_visc ! The time step for a part of the update due to viscosity [T ~> s]. @@ -240,6 +244,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & h_av(:,:,:) = 0; hp(:,:,:) = 0 up(:,:,:) = 0; upp(:,:,:) = 0 vp(:,:,:) = 0; vpp(:,:,:) = 0 + ueffA(:,:,:) = 0; veffA(:,:,:) = 0 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end) if (dyn_p_surf) then @@ -265,7 +270,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = u*h ! hp = h + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh, vh, dt*0.5, G, GV, US, CS%continuity_CSp, CS%OBC) + call continuity(u, v, h, hp, uh, vh, dt*0.5, G, GV, US, CS%continuity_CSp, CS%OBC, pbv) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -302,7 +307,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! CAu = -(f+zeta)/h_av vh + d/dx KE call cpu_clock_begin(id_clock_Cor) call CorAdCalc(u, v, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, GV, US, CS%CoriolisAdv_CSp) + G, GV, US, CS%CoriolisAdv_CSp, pbv) call cpu_clock_end(id_clock_Cor) ! PFu = d/dx M(h_av,T,S) @@ -355,7 +360,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = up * hp ! h_av = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, CS%continuity_CSp, CS%OBC) + call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, CS%continuity_CSp, CS%OBC, pbv) call cpu_clock_end(id_clock_continuity) call pass_var(h_av, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -368,7 +373,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! CAu = -(f+zeta(up))/h_av vh + d/dx KE(up) call cpu_clock_begin(id_clock_Cor) call CorAdCalc(up, vp, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, GV, US, CS%CoriolisAdv_CSp) + G, GV, US, CS%CoriolisAdv_CSp, pbv) call cpu_clock_end(id_clock_Cor) ! PFu = d/dx M(h_av,T,S) @@ -415,7 +420,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = upp * hp ! h = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, CS%continuity_CSp, CS%OBC) + call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, CS%continuity_CSp, CS%OBC, pbv) call cpu_clock_end(id_clock_continuity) call pass_var(h, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -430,6 +435,22 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & call disable_averaging(CS%diag) call enable_averages(dt, Time_local, CS%diag) +! Calculate effective areas and post data + if (CS%id_ueffA > 0) then + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_ueffA, ueffA, CS%diag) + endif + + if (CS%id_veffA > 0) then + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_veffA, veffA, CS%diag) + endif + + ! h_av = (h + hp)/2 do k=1,nz do j=js-2,je+2 ; do i=is-2,ie+2 @@ -446,7 +467,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! CAu = -(f+zeta(upp))/h_av vh + d/dx KE(upp) call cpu_clock_begin(id_clock_Cor) call CorAdCalc(upp, vpp, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, GV, US, CS%CoriolisAdv_CSp) + G, GV, US, CS%CoriolisAdv_CSp, pbv) call cpu_clock_end(id_clock_Cor) ! PFu = d/dx M(h_av,T,S) @@ -684,6 +705,12 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS 'Zonal Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & + 'Effective U Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + y_cell_method='sum', v_extensive = .true.) + CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & + 'Effective V Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + x_cell_method='sum', v_extensive = .true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index e6fec7f61e..a7d9abc856 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -48,7 +48,7 @@ module MOM_dynamics_unsplit_RK2 !* * !********+*********+*********+*********+*********+*********+*********+** -use MOM_variables, only : vertvisc_type, thermo_var_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_ptrs use MOM_variables, only : ocean_internal_state, accel_diag_ptrs, cont_diag_ptrs use MOM_forcing_type, only : mech_forcing use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum @@ -130,6 +130,7 @@ module MOM_dynamics_unsplit_RK2 !>@{ Diagnostic IDs integer :: id_uh = -1, id_vh = -1 + integer :: id_ueffA = -1, id_veffA = -1 integer :: id_PFu = -1, id_PFv = -1, id_CAu = -1, id_CAv = -1 !>@} @@ -188,7 +189,7 @@ module MOM_dynamics_unsplit_RK2 !> Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE) + VarMix, MEKE, pbv) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -234,10 +235,13 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, type(MEKE_type), pointer :: MEKE !< A pointer to a structure containing !! fields related to the Mesoscale !! Eddy Kinetic Energy. + type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av, hp real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocities [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocities [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffA ! Effective Area of U-Faces [H L ~> m2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffA ! Effective Area of V-Faces [H L ~> m2] real, dimension(:,:), pointer :: p_surf => NULL() real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s] real :: dt_visc ! The time step for a part of the update due to viscosity [T ~> s] @@ -250,6 +254,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, h_av(:,:,:) = 0; hp(:,:,:) = 0 up(:,:,:) = 0 vp(:,:,:) = 0 + ueffA(:,:,:) = 0; veffA(:,:,:) = 0 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end) if (dyn_p_surf) then @@ -281,7 +286,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call cpu_clock_begin(id_clock_continuity) ! This is a duplicate calculation of the last continuity from the previous step ! and could/should be optimized out. -AJA - call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, CS%continuity_CSp, CS%OBC) + call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, CS%continuity_CSp, CS%OBC, pbv) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -296,7 +301,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! CAu = -(f+zeta)/h_av vh + d/dx KE (function of u[n-1] and uh[n-1]) call cpu_clock_begin(id_clock_Cor) call CorAdCalc(u_in, v_in, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, GV, US, CS%CoriolisAdv_CSp) + G, GV, US, CS%CoriolisAdv_CSp, pbv) call cpu_clock_end(id_clock_Cor) ! PFu = d/dx M(h_av,T,S) (function of h[n-1/2]) @@ -350,7 +355,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! uh = up[n-1/2] * h[n-1/2] ! h_av = h + dt div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h_in, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC) + call continuity(up, vp, h_in, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -366,7 +371,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! CAu = -(f+zeta(up))/h_av vh + d/dx KE(up) (function of up[n-1/2], h[n-1/2]) call cpu_clock_begin(id_clock_Cor) call CorAdCalc(up, vp, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, & - G, GV, US, CS%CoriolisAdv_CSp) + G, GV, US, CS%CoriolisAdv_CSp, pbv) call cpu_clock_end(id_clock_Cor) if (associated(CS%OBC)) then call open_boundary_zero_normal_flow(CS%OBC, G, GV, CS%CAu, CS%CAv) @@ -405,7 +410,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! uh = up[n] * h[n] (up[n] might be extrapolated to damp GWs) ! h[n+1] = h[n] + dt div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h_in, h_in, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC) + call continuity(up, vp, h_in, h_in, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv) call cpu_clock_end(id_clock_continuity) call pass_var(h_in, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -448,6 +453,22 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, if (CS%id_uh > 0) call post_data(CS%id_uh, uh, CS%diag) if (CS%id_vh > 0) call post_data(CS%id_vh, vh, CS%diag) +! Calculate effective areas and post data + if (CS%id_ueffA > 0) then + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k)/up(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_ueffA, ueffA, CS%diag) + endif + + if (CS%id_veffA > 0) then + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k)/vp(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_veffA, veffA, CS%diag) + endif + + end subroutine step_MOM_dyn_unsplit_RK2 ! ============================================================================= @@ -645,6 +666,12 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag 'Zonal Pressure Force Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) + CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & + 'Effective U-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + y_cell_method='sum', v_extensive = .true.) + CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & + 'Effective V-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + x_cell_method='sum', v_extensive = .true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 7592dc8477..90482c6754 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -112,6 +112,16 @@ module MOM_grid IareaCv, & !< The masked inverse areas of v-grid cells [L-2 ~> m-2]. areaCv !< The areas of the v-grid cells [L2 ~> m2]. + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & + porous_DminU, & !< minimum topographic height of U-face [m] + porous_DmaxU, & !< maximum topographic height of U-face [m] + porous_DavgU !< average topographic height of U-face [m] + + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & + porous_DminV, & !< minimum topographic height of V-face [m] + porous_DmaxV, & !< maximum topographic height of V-face [m] + porous_DavgV !< average topographic height of V-face [m] + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim]. geoLatBu, & !< The geographic latitude at q points in degrees of latitude or m. @@ -574,6 +584,14 @@ subroutine allocate_metrics(G) ALLOC_(G%dx_Cv(isd:ied,JsdB:JedB)) ; G%dx_Cv(:,:) = 0.0 ALLOC_(G%dy_Cu(IsdB:IedB,jsd:jed)) ; G%dy_Cu(:,:) = 0.0 + ALLOC_(G%porous_DminU(IsdB:IedB,jsd:jed)); G%porous_DminU(:,:) = 0.0 + ALLOC_(G%porous_DmaxU(IsdB:IedB,jsd:jed)); G%porous_DmaxU(:,:) = 0.0 + ALLOC_(G%porous_DavgU(IsdB:IedB,jsd:jed)); G%porous_DavgU(:,:) = 0.0 + + ALLOC_(G%porous_DminV(isd:ied,JsdB:JedB)); G%porous_DminV(:,:) = 0.0 + ALLOC_(G%porous_DmaxV(isd:ied,JsdB:JedB)); G%porous_DmaxV(:,:) = 0.0 + ALLOC_(G%porous_DavgV(isd:ied,JsdB:JedB)); G%porous_DavgV(:,:) = 0.0 + ALLOC_(G%areaCu(IsdB:IedB,jsd:jed)) ; G%areaCu(:,:) = 0.0 ALLOC_(G%areaCv(isd:ied,JsdB:JedB)) ; G%areaCv(:,:) = 0.0 ALLOC_(G%IareaCu(IsdB:IedB,jsd:jed)) ; G%IareaCu(:,:) = 0.0 @@ -630,6 +648,9 @@ subroutine MOM_grid_end(G) DEALLOC_(G%dF_dx) ; DEALLOC_(G%dF_dy) DEALLOC_(G%sin_rot) ; DEALLOC_(G%cos_rot) + DEALLOC_(G%porous_DminU) ; DEALLOC_(G%porous_DmaxU) ; DEALLOC_(G%porous_DavgU) + DEALLOC_(G%porous_DminV) ; DEALLOC_(G%porous_DmaxV) ; DEALLOC_(G%porous_DavgV) + deallocate(G%gridLonT) ; deallocate(G%gridLatT) deallocate(G%gridLonB) ; deallocate(G%gridLatB) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 new file mode 100644 index 0000000000..4220f2c462 --- /dev/null +++ b/src/core/MOM_porous_barriers.F90 @@ -0,0 +1,166 @@ +!> Function for calculating curve fit for porous topography. +!written by sjd +module MOM_porous_barriers + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_error_handler, only : MOM_error, FATAL +use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs, porous_barrier_ptrs +use MOM_verticalGrid, only : verticalGrid_type +use MOM_interface_heights, only : find_eta + +implicit none ; private + +#include + +public porous_widths + +!> Calculates curve fit from D_min, D_max, D_avg +interface porous_widths + module procedure por_widths, calc_por_layer +end interface porous_widths + +contains + +subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) + !eta_bt, halo_size, eta_to_m not currently used + !variables needed to call find_eta + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: eta !< layer interface heights + !! [Z ~> m] or 1/eta_to_m m). + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic + !! variable that gives the "correct" free surface height (Boussinesq) or total water + !! column mass per unit area (non-Boussinesq). This is used to dilate the layer. + !! thicknesses when calculating interfaceheights [H ~> m or kg m-2]. + integer, optional, intent(in) :: halo_size !< width of halo points on + !! which to calculate eta. + + real, optional, intent(in) :: eta_to_m !< The conversion factor from + !! the units of eta to m; by default this is US%Z_to_m. + type(porous_barrier_ptrs), intent(inout) :: pbv !< porous barrier fractional cell metrics + + !local variables + integer ii, i, j, k, nk, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB + real w_layer, & ! fractional open width of layer interface [nondim] + A_layer, & ! integral of fractional open width from bottom to current layer[nondim] + A_layer_prev, & ! integral of fractional open width from bottom to previous layer [nondim] + eta_s, & ! layer height used for fit [Z ~> m] + eta_prev ! interface height of previous layer [Z ~> m] + isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed + IsdB = G%IsdB; IedB = G%IedB; JsdB = G%JsdB; JedB = G%JedB + + !eta is zero at surface and decreases downward + !all calculations are done in [m] + + nk = SZK_(G) + + !currently no treatment for using optional find_eta arguments if present + call find_eta(h, tv, G, GV, US, eta) + + do I=IsdB,IedB; do j=jsd,jed + if (G%porous_DavgU(I,j) < 0.) then + do K = nk+1,1,-1 + eta_s = max(US%Z_to_m*eta(I,j,K), US%Z_to_m*eta(I+1,j,K)) !take shallower layer height + !eta_s = 0.5 * (US%Z_to_m*eta(I,j,K) + US%Z_to_m*eta(I+1,j,K)) !take arithmetic mean + if (eta_s <= G%porous_DminU(I,j)) then + pbv%por_layer_widthU(I,j,K) = 0.0 + A_layer_prev = 0.0 + if (K < nk+1) then + pbv%por_face_areaU(I,j,k) = 0.0; endif + else + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j),& + eta_s, w_layer, A_layer) + pbv%por_layer_widthU(I,j,K) = w_layer + if (k <= nk) then + if ((eta_s - eta_prev) > 0.0) then + pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/& + (eta_s-eta_prev) + else + pbv%por_face_areaU(I,j,k) = 0.0; endif + endif + eta_prev = eta_s + A_layer_prev = A_layer + endif; enddo + endif; enddo; enddo + + do i=isd,ied; do J=JsdB,JedB + if (G%porous_DavgV(i,J) < 0.) then + do K = nk+1,1,-1 + eta_s = max(US%Z_to_m*eta(i,J,K), US%Z_to_m*eta(i,J+1,K)) !take shallower layer height + !eta_s = 0.5 * (US%Z_to_m*eta(i,J,K) + US%Z_to_m*eta(i,J+1,K)) !take arithmetic mean + if (eta_s <= G%porous_DminV(i,J)) then + pbv%por_layer_widthV(i,J,K) = 0.0 + A_layer_prev = 0.0 + if (K < nk+1) then + pbv%por_face_areaV(i,J,k) = 0.0; endif + else + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J),& + eta_s, w_layer, A_layer) + pbv%por_layer_widthV(i,J,K) = w_layer + if (k <= nk) then + if ((eta_s - eta_prev) > 0.0) then + pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/& + (eta_s-eta_prev) + else + pbv%por_face_areaU(I,j,k) = 0.0; endif + endif + eta_prev = eta_s + A_layer_prev = A_layer + endif; enddo + endif; enddo; enddo + +end subroutine por_widths + +subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) +!subroutine to calculate the profile fit for a layer + + real, intent(in) :: D_min !< minimum topographic height [m] + real, intent(in) :: D_max !< maximum topographic height [m] + real, intent(in) :: D_avg !< mean topographic height [m] + real, intent(in) :: eta_layer !< height of interface [m] + real, intent(out) :: w_layer !< frac. open interface width of current layer [nondim] + real, intent(out) :: A_layer !< frac. open face area of current layer [nondim] + !local variables + real m, a, & !convenience constant for fit [nondim] + zeta, & !normalized vertical coordinate [nondim] + psi, & !fractional width of layer between Dmin and Dmax [nondim] + psi_int !integral of psi from 0 to zeta + + !three parameter fit from Adcroft 2013 + m = (D_avg - D_min)/(D_max - D_min) + a = (1. - m)/m + + zeta = (eta_layer - D_min)/(D_max - D_min) + + if (eta_layer <= D_min) then + w_layer = 0.0 + A_layer = 0.0 + elseif (eta_layer >= D_max) then + w_layer = 1.0 + A_layer = eta_layer - D_avg + else + if (m < 0.5) then + psi = zeta**(1./a) + psi_int = (1.-m)*zeta**(1./(1.-m)) + elseif (m == 0.5) then + psi = zeta + psi_int = 0.5*zeta*zeta + else + psi = 1. - (1. - zeta)**a + psi_int = zeta - m + m*((1-zeta)**(1/m)) + endif + w_layer = psi + A_layer = (D_max - D_min)*psi_int + endif + + +end subroutine calc_por_layer + +end module MOM_porous_barriers diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index d3447f6590..e19df5b6c6 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -71,6 +71,10 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%dyCu(I,j) = dG%dyCu(I+ido,j+jdo) oG%dy_Cu(I,j) = dG%dy_Cu(I+ido,j+jdo) + oG%porous_DminU(I,j) = dG%porous_DminU(I+ido,j+jdo) + oG%porous_DmaxU(I,j) = dG%porous_DmaxU(I+ido,j+jdo) + oG%porous_DavgU(I,j) = dG%porous_DavgU(I+ido,j+jdo) + oG%mask2dCu(I,j) = dG%mask2dCu(I+ido,j+jdo) oG%areaCu(I,j) = dG%areaCu(I+ido,j+jdo) oG%IareaCu(I,j) = dG%IareaCu(I+ido,j+jdo) @@ -83,6 +87,10 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%dyCv(i,J) = dG%dyCv(i+ido,J+jdo) oG%dx_Cv(i,J) = dG%dx_Cv(i+ido,J+jdo) + oG%porous_DminV(i,J) = dG%porous_DminV(i+ido,J+jdo) + oG%porous_DmaxV(i,J) = dG%porous_DmaxV(i+ido,J+jdo) + oG%porous_DavgV(i,J) = dG%porous_DavgV(i+ido,J+jdo) + oG%mask2dCv(i,J) = dG%mask2dCv(i+ido,J+jdo) oG%areaCv(i,J) = dG%areaCv(i+ido,J+jdo) oG%IareaCv(i,J) = dG%IareaCv(i+ido,J+jdo) @@ -216,6 +224,10 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%dyCu(I,j) = oG%dyCu(I+ido,j+jdo) dG%dy_Cu(I,j) = oG%dy_Cu(I+ido,j+jdo) + dG%porous_DminU(I,j) = oG%porous_DminU(I+ido,j+jdo) + dG%porous_DmaxU(I,j) = oG%porous_DmaxU(I+ido,j+jdo) + dG%porous_DavgU(I,j) = oG%porous_DavgU(I+ido,j+jdo) + dG%mask2dCu(I,j) = oG%mask2dCu(I+ido,j+jdo) dG%areaCu(I,j) = oG%areaCu(I+ido,j+jdo) dG%IareaCu(I,j) = oG%IareaCu(I+ido,j+jdo) @@ -228,6 +240,10 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%dyCv(i,J) = oG%dyCv(i+ido,J+jdo) dG%dx_Cv(i,J) = oG%dx_Cv(i+ido,J+jdo) + dG%porous_DminV(i,J) = oG%porous_DminU(i+ido,J+jdo) + dG%porous_DmaxV(i,J) = oG%porous_DmaxU(i+ido,J+jdo) + dG%porous_DavgV(i,J) = oG%porous_DavgU(i+ido,J+jdo) + dG%mask2dCv(i,J) = oG%mask2dCv(i+ido,J+jdo) dG%areaCv(i,J) = oG%areaCv(i+ido,J+jdo) dG%IareaCv(i,J) = oG%IareaCv(i+ido,J+jdo) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 363f3eebfb..aff5a5438f 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -304,6 +304,16 @@ module MOM_variables type(group_pass_type) :: pass_FA_uv !< Structure for face area group halo updates end type BT_cont_type + +!> pointers to grids modifying cell metric at porous barriers +type, public :: porous_barrier_ptrs + real, pointer, dimension(:,:,:) :: por_face_areaU => NULL() !< fractional open area of U-faces [nondim] + real, pointer, dimension(:,:,:) :: por_face_areaV => NULL() !< fractional open area of V-faces [nondim] + real, pointer, dimension(:,:,:) :: por_layer_widthU => NULL() !< fractional open width of U-faces [nondim] + real, pointer, dimension(:,:,:) :: por_layer_widthV => NULL() !< fractional open width of V-faces [nondim] +end type porous_barrier_ptrs + + contains !> Allocates the fields for the surface (return) properties of diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index c7db67ee17..49d3dc01c7 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -109,6 +109,16 @@ module MOM_dyn_horgrid IareaCv, & !< The masked inverse areas of v-grid cells [L-2 ~> m-2]. areaCv !< The areas of the v-grid cells [L2 ~> m2]. + real, allocatable, dimension(:,:) :: & + porous_DminU, & !< minimum topographic height of U-face [m] + porous_DmaxU, & !< maximum topographic height of U-face [m] + porous_DavgU !< average topographic height of U-face [m] + + real, allocatable, dimension(:,:) :: & + porous_DminV, & !< minimum topographic height of V-face [m] + porous_DmaxV, & !< maximum topographic height of V-face [m] + porous_DavgV !< average topographic height of V-face [m] + real, allocatable, dimension(:,:) :: & mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim]. geoLatBu, & !< The geographic latitude at q points [degrees of latitude] or [m]. @@ -256,6 +266,15 @@ subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel) allocate(G%IareaCu(IsdB:IedB,jsd:jed), source=0.0) allocate(G%IareaCv(isd:ied,JsdB:JedB), source=0.0) + allocate(G%porous_DminU(IsdB:IedB,jsd:jed), source=0.0) + allocate(G%porous_DmaxU(IsdB:IedB,jsd:jed), source=0.0) + allocate(G%porous_DavgU(IsdB:IedB,jsd:jed), source=0.0) + + allocate(G%porous_DminV(isd:ied,JsdB:JedB), source=0.0) + allocate(G%porous_DmaxV(isd:ied,JsdB:JedB), source=0.0) + allocate(G%porous_DavgV(isd:ied,JsdB:JedB), source=0.0) + + allocate(G%bathyT(isd:ied, jsd:jed), source=0.0) allocate(G%CoriolisBu(IsdB:IedB, JsdB:JedB), source=0.0) allocate(G%dF_dx(isd:ied, jsd:jed), source=0.0) @@ -489,6 +508,9 @@ subroutine destroy_dyn_horgrid(G) deallocate(G%dx_Cv) ; deallocate(G%dy_Cu) + deallocate(G%porous_DminU) ; deallocate(G%porous_DmaxU) ; deallocate(G%porous_DavgU) + deallocate(G%porous_DminV) ; deallocate(G%porous_DmaxV) ; deallocate(G%porous_DavgV) + deallocate(G%bathyT) ; deallocate(G%CoriolisBu) deallocate(G%dF_dx) ; deallocate(G%dF_dy) deallocate(G%sin_rot) ; deallocate(G%cos_rot) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index c252e296a5..e9531bb4e2 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -850,6 +850,10 @@ subroutine reset_face_lengths_list(G, param_file, US) integer, allocatable, dimension(:) :: & u_line_no, v_line_no, & ! The line numbers in lines of u- and v-face lines u_line_used, v_line_used ! The number of times each u- and v-line is used. + real, allocatable, dimension(:) :: & + Dmin_u, Dmax_u, Davg_u ! Porous barrier monomial fit params [m] + real, allocatable, dimension(:) :: & + Dmin_v, Dmax_v, Davg_v real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] real :: lat, lon ! The latitude and longitude of a point. @@ -865,6 +869,9 @@ subroutine reset_face_lengths_list(G, param_file, US) integer :: ios, iounit, isu, isv integer :: last, num_lines, nl_read, ln, npt, u_pt, v_pt integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB + integer :: isu_por, isv_por + logical :: found_u_por, found_v_por + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -929,6 +936,14 @@ subroutine reset_face_lengths_list(G, param_file, US) allocate(v_line_used(num_lines), source=0) allocate(v_line_no(num_lines), source=0) + allocate(Dmin_u(num_lines)) ; Dmin_u(:) = 0.0 + allocate(Dmax_u(num_lines)) ; Dmax_u(:) = 0.0 + allocate(Davg_u(num_lines)) ; Davg_u(:) = 0.0 + + allocate(Dmin_v(num_lines)) ; Dmin_v(:) = 0.0 + allocate(Dmax_v(num_lines)) ; Dmax_v(:) = 0.0 + allocate(Davg_v(num_lines)) ; Davg_v(:) = 0.0 + ! Actually read the lines. if (is_root_pe()) then call read_face_length_list(iounit, filename, nl_read, lines) @@ -946,13 +961,21 @@ subroutine reset_face_lengths_list(G, param_file, US) line = lines(ln) ! Detect keywords found_u = .false.; found_v = .false. + found_u_por = .false.; found_v_por = .false. isu = index(uppercase(line), "U_WIDTH" ); if (isu > 0) found_u = .true. isv = index(uppercase(line), "V_WIDTH" ); if (isv > 0) found_v = .true. + isu_por = index(uppercase(line), "U_WIDTH_POR" ); if (isu_por > 0) found_u_por = .true. + isv_por = index(uppercase(line), "V_WIDTH_POR" ); if (isv_por > 0) found_v_por = .true. ! Store and check the relevant values. if (found_u) then u_pt = u_pt + 1 - read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt) + if (found_u_por .eqv. .false.) then + read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt) + elseif (found_u_por) then + read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), & + Dmin_u(u_pt), Dmax_u(u_pt), Davg_u(u_pt) + endif u_line_no(u_pt) = ln if (is_root_PE()) then if (check_360) then @@ -977,10 +1000,19 @@ subroutine reset_face_lengths_list(G, param_file, US) call MOM_error(WARNING, "reset_face_lengths_list : Negative "//& "u-width found when reading line "//trim(line)//" from file "//& trim(filename)) + if (Dmin_u(u_pt) > Dmax_u(u_pt)) & + call MOM_error(WARNING, "reset_face_lengths_list : Out-of-order "//& + "topographical min/max found when reading line "//trim(line)//" from file "//& + trim(filename)) endif elseif (found_v) then v_pt = v_pt + 1 - read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt) + if (found_v_por .eqv. .false.) then + read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt) + elseif (found_v_por) then + read(line(isv+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), & + Dmin_v(v_pt), Dmax_v(v_pt), Davg_v(v_pt) + endif v_line_no(v_pt) = ln if (is_root_PE()) then if (check_360) then @@ -1005,6 +1037,10 @@ subroutine reset_face_lengths_list(G, param_file, US) call MOM_error(WARNING, "reset_face_lengths_list : Negative "//& "v-width found when reading line "//trim(line)//" from file "//& trim(filename)) + if (Dmin_v(v_pt) > Dmax_v(v_pt)) & + call MOM_error(WARNING, "reset_face_lengths_list : Out-of-order "//& + "topographical min/max found when reading line "//trim(line)//" from file "//& + trim(filename)) endif endif enddo @@ -1023,6 +1059,10 @@ subroutine reset_face_lengths_list(G, param_file, US) ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then G%dy_Cu(I,j) = G%mask2dCu(I,j) * m_to_L*min(L_to_m*G%dyCu(I,j), max(u_width(npt), 0.0)) + G%porous_DminU(I,j) = Dmin_u(npt) + G%porous_DmaxU(I,j) = Dmax_u(npt) + G%porous_DavgU(I,j) = Davg_u(npt) + if (j>=G%jsc .and. j<=G%jec .and. I>=G%isc .and. I<=G%iec) then ! Limit messages/checking to compute domain if ( G%mask2dCu(I,j) == 0.0 ) then write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",& @@ -1032,6 +1072,9 @@ subroutine reset_face_lengths_list(G, param_file, US) write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",& u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") to ",L_to_m*G%dy_Cu(I,j),"m" + write(stdout,'(A,3F8.2,A)') & + "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",G%porous_DminU(I,j),& + G%porous_DmaxU(I,j), G%porous_DavgU(I,j),")m" endif endif endif @@ -1053,6 +1096,9 @@ subroutine reset_face_lengths_list(G, param_file, US) ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. & ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then G%dx_Cv(i,J) = G%mask2dCv(i,J) * m_to_L*min(L_to_m*G%dxCv(i,J), max(v_width(npt), 0.0)) + G%porous_DminV(i,J) = Dmin_v(npt) + G%porous_DmaxV(i,J) = Dmax_v(npt) + G%porous_DavgV(i,J) = Davg_v(npt) if (i>=G%isc .and. i<=G%iec .and. J>=G%jsc .and. J<=G%jec) then ! Limit messages/checking to compute domain if ( G%mask2dCv(i,J) == 0.0 ) then write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",& @@ -1062,6 +1108,9 @@ subroutine reset_face_lengths_list(G, param_file, US) write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",& v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") to ",L_to_m*G%dx_Cv(I,j),"m" + write(stdout,'(A,3F8.2,A)') & + "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",G%porous_DminV(i,J),& + G%porous_DmaxV(i,J), G%porous_DavgV(i,J),")m" endif endif endif @@ -1097,6 +1146,8 @@ subroutine reset_face_lengths_list(G, param_file, US) deallocate(u_line_used, v_line_used, u_line_no, v_line_no) deallocate(u_lat) ; deallocate(u_lon) ; deallocate(u_width) deallocate(v_lat) ; deallocate(v_lon) ; deallocate(v_width) + deallocate(Dmin_u) ; deallocate(Dmax_u) ; deallocate(Davg_u) + deallocate(Dmin_v) ; deallocate(Dmax_v) ; deallocate(Davg_v) endif call callTree_leave(trim(mdl)//'()') diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 1cf3b5ddc9..60d9d49d7c 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -23,7 +23,7 @@ module MOM_set_visc use MOM_restart, only : register_restart_field_as_obsolete use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs, vertvisc_type +use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_DIRECTION_E @@ -115,7 +115,7 @@ module MOM_set_visc contains !> Calculates the thickness of the bottom boundary layer and the viscosity within that layer. -subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) +subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -132,6 +132,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) !! related fields. type(set_visc_CS), pointer :: CS !< The control structure returned by a previous !! call to set_visc_init. + type(porous_barrier_ptrs),intent(in) :: pbv !< porous barrier fractional cell metrics ! Local variables real, dimension(SZIB_(G)) :: & @@ -373,7 +374,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) !$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,nz,nkmb, & !$OMP nkml,Isq,Ieq,Jsq,Jeq,h_neglect,Rho0x400_G,C2pi_3, & !$OMP U_bg_sq,cdrag_sqrt_Z,cdrag_sqrt,K2,use_BBL_EOS, & - !$OMP OBC,maxitt,D_u,D_v,mask_u,mask_v) & + !$OMP OBC,maxitt,D_u,D_v,mask_u,mask_v, pbv) & !$OMP firstprivate(Vol_quit) do j=Jsq,Jeq ; do m=1,2 @@ -903,6 +904,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) endif ! end of a<0 cases. endif + !modify L(K) for porous barrier parameterization + if (m==1) then ; L(K) = L(K)*pbv%por_layer_widthU(I,j,K) + else ; L(K) = L(K)*pbv%por_layer_widthV(i,J,K); endif + ! Determine the drag contributing to the bottom boundary layer ! and the Raleigh drag that acts on each layer. if (L(K) > L(K+1)) then @@ -913,8 +918,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) BBL_frac = 0.0 endif - if (m==1) then ; Cell_width = G%dy_Cu(I,j) - else ; Cell_width = G%dx_Cv(i,J) ; endif + if (m==1) then ; Cell_width = G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k) + else ; Cell_width = G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k) ; endif gam = 1.0 - L(K+1)/L(K) Rayleigh = US%L_to_Z * CS%cdrag * (L(K)-L(K+1)) * (1.0-BBL_frac) * & (12.0*CS%c_Smag*h_vel_pos) / (12.0*CS%c_Smag*h_vel_pos + & From bc826090c38f372e459f7b35d122baf81cd16218 Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Fri, 19 Nov 2021 14:57:00 -0500 Subject: [PATCH 2/8] small fixes to porous topo code * Style Fixes * Store depths for porous curve fitting as [Z ~> m] and account for Z_ref --- src/core/MOM_dynamics_split_RK2.F90 | 3 +- src/core/MOM_dynamics_unsplit.F90 | 3 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 3 +- src/core/MOM_grid.F90 | 12 +- src/core/MOM_porous_barriers.F90 | 138 +++++++++--------- src/core/MOM_transcribe_grid.F90 | 24 +-- src/framework/MOM_dyn_horgrid.F90 | 12 +- .../MOM_shared_initialization.F90 | 22 +-- 8 files changed, 114 insertions(+), 103 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index fa796c7fe8..e3f30417ef 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -398,7 +398,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s do j=G%jsdB,G%jedB ; do i=G%isd,G%ied ; vp(i,j,k) = 0.0 ; enddo ; enddo do j=G%jsd,G%jed ; do i=G%isd,G%ied ; hp(i,j,k) = h(i,j,k) ; enddo ; enddo enddo - ueffA(:,:,:) = 0; veffA(:,:,:) = 0 + if (CS%id_ueffA > 0) ueffA(:,:,:) = 0 + if (CS%id_veffA > 0) veffA(:,:,:) = 0 ! Update CFL truncation value as function of time call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 7cfc9d649c..e52085c2d5 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -244,7 +244,8 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & h_av(:,:,:) = 0; hp(:,:,:) = 0 up(:,:,:) = 0; upp(:,:,:) = 0 vp(:,:,:) = 0; vpp(:,:,:) = 0 - ueffA(:,:,:) = 0; veffA(:,:,:) = 0 + if (CS%id_ueffA > 0) ueffA(:,:,:) = 0 + if (CS%id_veffA > 0) veffA(:,:,:) = 0 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end) if (dyn_p_surf) then diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index a7d9abc856..f8112fe4cd 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -254,7 +254,8 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, h_av(:,:,:) = 0; hp(:,:,:) = 0 up(:,:,:) = 0 vp(:,:,:) = 0 - ueffA(:,:,:) = 0; veffA(:,:,:) = 0 + if (CS%id_ueffA > 0) ueffA(:,:,:) = 0 + if (CS%id_veffA > 0) veffA(:,:,:) = 0 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end) if (dyn_p_surf) then diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 90482c6754..b7ef04b76f 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -113,14 +113,14 @@ module MOM_grid areaCv !< The areas of the v-grid cells [L2 ~> m2]. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & - porous_DminU, & !< minimum topographic height of U-face [m] - porous_DmaxU, & !< maximum topographic height of U-face [m] - porous_DavgU !< average topographic height of U-face [m] + porous_DminU, & !< minimum topographic height of U-face [Z ~> m] + porous_DmaxU, & !< maximum topographic height of U-face [Z ~> m] + porous_DavgU !< average topographic height of U-face [Z ~> m] real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & - porous_DminV, & !< minimum topographic height of V-face [m] - porous_DmaxV, & !< maximum topographic height of V-face [m] - porous_DavgV !< average topographic height of V-face [m] + porous_DminV, & !< minimum topographic height of V-face [Z ~> m] + porous_DmaxV, & !< maximum topographic height of V-face [Z ~> m] + porous_DavgV !< average topographic height of V-face [Z ~> m] real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim]. diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 4220f2c462..1230b47cfb 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -1,4 +1,4 @@ -!> Function for calculating curve fit for porous topography. +!> Module for calculating curve fit for porous topography. !written by sjd module MOM_porous_barriers @@ -24,6 +24,7 @@ module MOM_porous_barriers contains +!> subroutine to assign cell face areas and layer widths for porous topography subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) !eta_bt, halo_size, eta_to_m not currently used !variables needed to call find_eta @@ -57,69 +58,74 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) IsdB = G%IsdB; IedB = G%IedB; JsdB = G%JsdB; JedB = G%JedB !eta is zero at surface and decreases downward - !all calculations are done in [m] nk = SZK_(G) !currently no treatment for using optional find_eta arguments if present call find_eta(h, tv, G, GV, US, eta) - do I=IsdB,IedB; do j=jsd,jed - if (G%porous_DavgU(I,j) < 0.) then - do K = nk+1,1,-1 - eta_s = max(US%Z_to_m*eta(I,j,K), US%Z_to_m*eta(I+1,j,K)) !take shallower layer height - !eta_s = 0.5 * (US%Z_to_m*eta(I,j,K) + US%Z_to_m*eta(I+1,j,K)) !take arithmetic mean - if (eta_s <= G%porous_DminU(I,j)) then - pbv%por_layer_widthU(I,j,K) = 0.0 - A_layer_prev = 0.0 - if (K < nk+1) then - pbv%por_face_areaU(I,j,k) = 0.0; endif - else - call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j),& - eta_s, w_layer, A_layer) - pbv%por_layer_widthU(I,j,K) = w_layer - if (k <= nk) then - if ((eta_s - eta_prev) > 0.0) then - pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/& + do j=jsd,jed; do I=IsdB,IedB + if (G%porous_DavgU(I,j) < 0.) then + do K = nk+1,1,-1 + eta_s = max(US%Z_to_m*eta(I,j,K), US%Z_to_m*eta(I+1,j,K)) !take shallower layer height + !eta_s = 0.5 * (US%Z_to_m*eta(I,j,K) + US%Z_to_m*eta(I+1,j,K)) !take arithmetic mean + if (eta_s <= G%porous_DminU(I,j)) then + pbv%por_layer_widthU(I,j,K) = 0.0 + A_layer_prev = 0.0 + if (K < nk+1) then + pbv%por_face_areaU(I,j,k) = 0.0; endif + else + call calc_por_layer(US%Z_to_m*(G%porous_DminU(I,j)-G%Z_ref), & + US%Z_to_m*(G%porous_DmaxU(I,j)-G%Z_ref), & + US%Z_to_m*(G%porous_DavgU(I,j)-G%Z_ref), eta_s, w_layer, A_layer) + pbv%por_layer_widthU(I,j,K) = w_layer + if (k <= nk) then + if ((eta_s - eta_prev) > 0.0) then + pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/& (eta_s-eta_prev) - else - pbv%por_face_areaU(I,j,k) = 0.0; endif - endif - eta_prev = eta_s - A_layer_prev = A_layer - endif; enddo - endif; enddo; enddo + else + pbv%por_face_areaU(I,j,k) = 0.0; endif + endif + eta_prev = eta_s + A_layer_prev = A_layer + endif + enddo + endif + enddo; enddo do i=isd,ied; do J=JsdB,JedB - if (G%porous_DavgV(i,J) < 0.) then - do K = nk+1,1,-1 - eta_s = max(US%Z_to_m*eta(i,J,K), US%Z_to_m*eta(i,J+1,K)) !take shallower layer height - !eta_s = 0.5 * (US%Z_to_m*eta(i,J,K) + US%Z_to_m*eta(i,J+1,K)) !take arithmetic mean - if (eta_s <= G%porous_DminV(i,J)) then - pbv%por_layer_widthV(i,J,K) = 0.0 - A_layer_prev = 0.0 - if (K < nk+1) then - pbv%por_face_areaV(i,J,k) = 0.0; endif - else - call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J),& - eta_s, w_layer, A_layer) - pbv%por_layer_widthV(i,J,K) = w_layer - if (k <= nk) then - if ((eta_s - eta_prev) > 0.0) then - pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/& + if (G%porous_DavgV(i,J) < 0.) then + do K = nk+1,1,-1 + eta_s = max(US%Z_to_m*eta(i,J,K), US%Z_to_m*eta(i,J+1,K)) !take shallower layer height + !eta_s = 0.5 * (US%Z_to_m*eta(i,J,K) + US%Z_to_m*eta(i,J+1,K)) !take arithmetic mean + if (eta_s <= G%porous_DminV(i,J)) then + pbv%por_layer_widthV(i,J,K) = 0.0 + A_layer_prev = 0.0 + if (K < nk+1) then + pbv%por_face_areaV(i,J,k) = 0.0; endif + else + call calc_por_layer(US%Z_to_m*(G%porous_DminV(i,J)-G%Z_ref), & + US%Z_to_m*(G%porous_DmaxV(i,J)-G%Z_ref), & + US%Z_to_m*(G%porous_DavgV(i,J)-G%Z_ref), eta_s, w_layer, A_layer) + pbv%por_layer_widthV(i,J,K) = w_layer + if (k <= nk) then + if ((eta_s - eta_prev) > 0.0) then + pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/& (eta_s-eta_prev) - else - pbv%por_face_areaU(I,j,k) = 0.0; endif - endif - eta_prev = eta_s - A_layer_prev = A_layer - endif; enddo - endif; enddo; enddo + else + pbv%por_face_areaU(I,j,k) = 0.0; endif + endif + eta_prev = eta_s + A_layer_prev = A_layer + endif + enddo + endif + enddo; enddo end subroutine por_widths +!> subroutine to calculate the profile fit for a single layer in a column subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) -!subroutine to calculate the profile fit for a layer real, intent(in) :: D_min !< minimum topographic height [m] real, intent(in) :: D_max !< maximum topographic height [m] @@ -140,24 +146,24 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) zeta = (eta_layer - D_min)/(D_max - D_min) if (eta_layer <= D_min) then - w_layer = 0.0 - A_layer = 0.0 + w_layer = 0.0 + A_layer = 0.0 elseif (eta_layer >= D_max) then - w_layer = 1.0 - A_layer = eta_layer - D_avg + w_layer = 1.0 + A_layer = eta_layer - D_avg else - if (m < 0.5) then - psi = zeta**(1./a) - psi_int = (1.-m)*zeta**(1./(1.-m)) - elseif (m == 0.5) then - psi = zeta - psi_int = 0.5*zeta*zeta - else - psi = 1. - (1. - zeta)**a - psi_int = zeta - m + m*((1-zeta)**(1/m)) - endif - w_layer = psi - A_layer = (D_max - D_min)*psi_int + if (m < 0.5) then + psi = zeta**(1./a) + psi_int = (1.-m)*zeta**(1./(1.-m)) + elseif (m == 0.5) then + psi = zeta + psi_int = 0.5*zeta*zeta + else + psi = 1. - (1. - zeta)**a + psi_int = zeta - m + m*((1-zeta)**(1/m)) + endif + w_layer = psi + A_layer = (D_max - D_min)*psi_int endif diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index e19df5b6c6..0e6d681896 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -71,9 +71,9 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%dyCu(I,j) = dG%dyCu(I+ido,j+jdo) oG%dy_Cu(I,j) = dG%dy_Cu(I+ido,j+jdo) - oG%porous_DminU(I,j) = dG%porous_DminU(I+ido,j+jdo) - oG%porous_DmaxU(I,j) = dG%porous_DmaxU(I+ido,j+jdo) - oG%porous_DavgU(I,j) = dG%porous_DavgU(I+ido,j+jdo) + oG%porous_DminU(I,j) = dG%porous_DminU(I+ido,j+jdo) - oG%Z_ref + oG%porous_DmaxU(I,j) = dG%porous_DmaxU(I+ido,j+jdo) - oG%Z_ref + oG%porous_DavgU(I,j) = dG%porous_DavgU(I+ido,j+jdo) - oG%Z_ref oG%mask2dCu(I,j) = dG%mask2dCu(I+ido,j+jdo) oG%areaCu(I,j) = dG%areaCu(I+ido,j+jdo) @@ -87,9 +87,9 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%dyCv(i,J) = dG%dyCv(i+ido,J+jdo) oG%dx_Cv(i,J) = dG%dx_Cv(i+ido,J+jdo) - oG%porous_DminV(i,J) = dG%porous_DminV(i+ido,J+jdo) - oG%porous_DmaxV(i,J) = dG%porous_DmaxV(i+ido,J+jdo) - oG%porous_DavgV(i,J) = dG%porous_DavgV(i+ido,J+jdo) + oG%porous_DminV(i,J) = dG%porous_DminV(i+ido,J+jdo) - oG%Z_ref + oG%porous_DmaxV(i,J) = dG%porous_DmaxV(i+ido,J+jdo) - oG%Z_ref + oG%porous_DavgV(i,J) = dG%porous_DavgV(i+ido,J+jdo) - oG%Z_ref oG%mask2dCv(i,J) = dG%mask2dCv(i+ido,J+jdo) oG%areaCv(i,J) = dG%areaCv(i+ido,J+jdo) @@ -224,9 +224,9 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%dyCu(I,j) = oG%dyCu(I+ido,j+jdo) dG%dy_Cu(I,j) = oG%dy_Cu(I+ido,j+jdo) - dG%porous_DminU(I,j) = oG%porous_DminU(I+ido,j+jdo) - dG%porous_DmaxU(I,j) = oG%porous_DmaxU(I+ido,j+jdo) - dG%porous_DavgU(I,j) = oG%porous_DavgU(I+ido,j+jdo) + dG%porous_DminU(I,j) = oG%porous_DminU(I+ido,j+jdo) + oG%Z_ref + dG%porous_DmaxU(I,j) = oG%porous_DmaxU(I+ido,j+jdo) + oG%Z_ref + dG%porous_DavgU(I,j) = oG%porous_DavgU(I+ido,j+jdo) + oG%Z_ref dG%mask2dCu(I,j) = oG%mask2dCu(I+ido,j+jdo) dG%areaCu(I,j) = oG%areaCu(I+ido,j+jdo) @@ -240,9 +240,9 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%dyCv(i,J) = oG%dyCv(i+ido,J+jdo) dG%dx_Cv(i,J) = oG%dx_Cv(i+ido,J+jdo) - dG%porous_DminV(i,J) = oG%porous_DminU(i+ido,J+jdo) - dG%porous_DmaxV(i,J) = oG%porous_DmaxU(i+ido,J+jdo) - dG%porous_DavgV(i,J) = oG%porous_DavgU(i+ido,J+jdo) + dG%porous_DminV(i,J) = oG%porous_DminU(i+ido,J+jdo) + oG%Z_ref + dG%porous_DmaxV(i,J) = oG%porous_DmaxU(i+ido,J+jdo) + oG%Z_ref + dG%porous_DavgV(i,J) = oG%porous_DavgU(i+ido,J+jdo) + oG%Z_ref dG%mask2dCv(i,J) = oG%mask2dCv(i+ido,J+jdo) dG%areaCv(i,J) = oG%areaCv(i+ido,J+jdo) diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 49d3dc01c7..efa3b02b2b 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -110,14 +110,14 @@ module MOM_dyn_horgrid areaCv !< The areas of the v-grid cells [L2 ~> m2]. real, allocatable, dimension(:,:) :: & - porous_DminU, & !< minimum topographic height of U-face [m] - porous_DmaxU, & !< maximum topographic height of U-face [m] - porous_DavgU !< average topographic height of U-face [m] + porous_DminU, & !< minimum topographic height of U-face [Z ~> m] + porous_DmaxU, & !< maximum topographic height of U-face [Z ~> m] + porous_DavgU !< average topographic height of U-face [Z ~> m] real, allocatable, dimension(:,:) :: & - porous_DminV, & !< minimum topographic height of V-face [m] - porous_DmaxV, & !< maximum topographic height of V-face [m] - porous_DavgV !< average topographic height of V-face [m] + porous_DminV, & !< minimum topographic height of V-face [Z ~> m] + porous_DmaxV, & !< maximum topographic height of V-face [Z ~> m] + porous_DavgV !< average topographic height of V-face [Z ~> m] real, allocatable, dimension(:,:) :: & mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim]. diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index e9531bb4e2..8f1a309d5f 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -856,6 +856,7 @@ subroutine reset_face_lengths_list(G, param_file, US) Dmin_v, Dmax_v, Davg_v real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + real :: m_to_Z ! A unit conversion factor [L ~> m] real :: lat, lon ! The latitude and longitude of a point. real :: len_lon ! The periodic range of longitudes, usually 360 degrees. real :: len_lat ! The range of latitudes, usually 180 degrees. @@ -878,6 +879,7 @@ subroutine reset_face_lengths_list(G, param_file, US) call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m + m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z call get_param(param_file, mdl, "CHANNEL_LIST_FILE", chan_file, & "The file from which the list of narrowed channels is read.", & @@ -971,9 +973,9 @@ subroutine reset_face_lengths_list(G, param_file, US) if (found_u) then u_pt = u_pt + 1 if (found_u_por .eqv. .false.) then - read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt) + read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt) elseif (found_u_por) then - read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), & + read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), & Dmin_u(u_pt), Dmax_u(u_pt), Davg_u(u_pt) endif u_line_no(u_pt) = ln @@ -1008,9 +1010,9 @@ subroutine reset_face_lengths_list(G, param_file, US) elseif (found_v) then v_pt = v_pt + 1 if (found_v_por .eqv. .false.) then - read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt) + read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt) elseif (found_v_por) then - read(line(isv+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), & + read(line(isv+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), & Dmin_v(v_pt), Dmax_v(v_pt), Davg_v(v_pt) endif v_line_no(v_pt) = ln @@ -1059,9 +1061,9 @@ subroutine reset_face_lengths_list(G, param_file, US) ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then G%dy_Cu(I,j) = G%mask2dCu(I,j) * m_to_L*min(L_to_m*G%dyCu(I,j), max(u_width(npt), 0.0)) - G%porous_DminU(I,j) = Dmin_u(npt) - G%porous_DmaxU(I,j) = Dmax_u(npt) - G%porous_DavgU(I,j) = Davg_u(npt) + G%porous_DminU(I,j) = m_to_Z*Dmin_u(npt) + G%porous_DmaxU(I,j) = m_to_Z*Dmax_u(npt) + G%porous_DavgU(I,j) = m_to_Z*Davg_u(npt) if (j>=G%jsc .and. j<=G%jec .and. I>=G%isc .and. I<=G%iec) then ! Limit messages/checking to compute domain if ( G%mask2dCu(I,j) == 0.0 ) then @@ -1096,9 +1098,9 @@ subroutine reset_face_lengths_list(G, param_file, US) ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. & ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then G%dx_Cv(i,J) = G%mask2dCv(i,J) * m_to_L*min(L_to_m*G%dxCv(i,J), max(v_width(npt), 0.0)) - G%porous_DminV(i,J) = Dmin_v(npt) - G%porous_DmaxV(i,J) = Dmax_v(npt) - G%porous_DavgV(i,J) = Davg_v(npt) + G%porous_DminV(i,J) = m_to_Z*Dmin_v(npt) + G%porous_DmaxV(i,J) = m_to_Z*Dmax_v(npt) + G%porous_DavgV(i,J) = m_to_Z*Davg_v(npt) if (i>=G%isc .and. i<=G%iec .and. J>=G%jsc .and. J<=G%jec) then ! Limit messages/checking to compute domain if ( G%mask2dCv(i,J) == 0.0 ) then write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",& From e717e46ee6593b4460db68b56938d62bd4ec05ec Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Wed, 24 Nov 2021 14:57:34 -0500 Subject: [PATCH 3/8] preserve dimensions in porous topo interpolation * changed porous topo code to preserve dimensions * readded rotations that got lost in merge --- src/core/MOM_porous_barriers.F90 | 67 ++++++++++++++++--------------- src/framework/MOM_dyn_horgrid.F90 | 8 ++++ 2 files changed, 43 insertions(+), 32 deletions(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 1230b47cfb..34af04a7e4 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -50,8 +50,8 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) !local variables integer ii, i, j, k, nk, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB real w_layer, & ! fractional open width of layer interface [nondim] - A_layer, & ! integral of fractional open width from bottom to current layer[nondim] - A_layer_prev, & ! integral of fractional open width from bottom to previous layer [nondim] + A_layer, & ! integral of fractional open width from bottom to current layer[Z ~> m] + A_layer_prev, & ! integral of fractional open width from bottom to previous layer [Z ~> m] eta_s, & ! layer height used for fit [Z ~> m] eta_prev ! interface height of previous layer [Z ~> m] isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed @@ -67,22 +67,20 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) do j=jsd,jed; do I=IsdB,IedB if (G%porous_DavgU(I,j) < 0.) then do K = nk+1,1,-1 - eta_s = max(US%Z_to_m*eta(I,j,K), US%Z_to_m*eta(I+1,j,K)) !take shallower layer height - !eta_s = 0.5 * (US%Z_to_m*eta(I,j,K) + US%Z_to_m*eta(I+1,j,K)) !take arithmetic mean - if (eta_s <= G%porous_DminU(I,j)) then + eta_s = max(eta(I,j,K), eta(I+1,j,K)) !take shallower layer height + if (US%Z_to_m*eta_s <= (US%Z_to_m*G%porous_DminU(I,j))) then pbv%por_layer_widthU(I,j,K) = 0.0 A_layer_prev = 0.0 if (K < nk+1) then pbv%por_face_areaU(I,j,k) = 0.0; endif else - call calc_por_layer(US%Z_to_m*(G%porous_DminU(I,j)-G%Z_ref), & - US%Z_to_m*(G%porous_DmaxU(I,j)-G%Z_ref), & - US%Z_to_m*(G%porous_DavgU(I,j)-G%Z_ref), eta_s, w_layer, A_layer) + call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), & + G%porous_DavgU(I,j), eta_s, w_layer, A_layer, G%Z_ref, US%Z_to_m) pbv%por_layer_widthU(I,j,K) = w_layer if (k <= nk) then if ((eta_s - eta_prev) > 0.0) then pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/& - (eta_s-eta_prev) + (US%Z_to_m*(eta_s-eta_prev)) else pbv%por_face_areaU(I,j,k) = 0.0; endif endif @@ -96,22 +94,20 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) do i=isd,ied; do J=JsdB,JedB if (G%porous_DavgV(i,J) < 0.) then do K = nk+1,1,-1 - eta_s = max(US%Z_to_m*eta(i,J,K), US%Z_to_m*eta(i,J+1,K)) !take shallower layer height - !eta_s = 0.5 * (US%Z_to_m*eta(i,J,K) + US%Z_to_m*eta(i,J+1,K)) !take arithmetic mean - if (eta_s <= G%porous_DminV(i,J)) then + eta_s = max(eta(i,J,K), eta(i,J+1,K)) !take shallower layer height + if (US%Z_to_m*eta_s <= US%Z_to_m*G%porous_DminV(i,J)) then pbv%por_layer_widthV(i,J,K) = 0.0 A_layer_prev = 0.0 if (K < nk+1) then pbv%por_face_areaV(i,J,k) = 0.0; endif else - call calc_por_layer(US%Z_to_m*(G%porous_DminV(i,J)-G%Z_ref), & - US%Z_to_m*(G%porous_DmaxV(i,J)-G%Z_ref), & - US%Z_to_m*(G%porous_DavgV(i,J)-G%Z_ref), eta_s, w_layer, A_layer) + call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), & + G%porous_DavgV(i,J), eta_s, w_layer, A_layer, G%Z_ref, US%Z_to_m) pbv%por_layer_widthV(i,J,K) = w_layer if (k <= nk) then if ((eta_s - eta_prev) > 0.0) then pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/& - (eta_s-eta_prev) + (US%Z_to_m*(eta_s-eta_prev)) else pbv%por_face_areaU(I,j,k) = 0.0; endif endif @@ -125,32 +121,39 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) end subroutine por_widths !> subroutine to calculate the profile fit for a single layer in a column -subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) +subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer, Z_ref, Z_to_m) - real, intent(in) :: D_min !< minimum topographic height [m] - real, intent(in) :: D_max !< maximum topographic height [m] - real, intent(in) :: D_avg !< mean topographic height [m] - real, intent(in) :: eta_layer !< height of interface [m] + real, intent(in) :: D_min !< minimum topographic height [Z ~> m] + real, intent(in) :: D_max !< maximum topographic height [Z ~> m] + real, intent(in) :: D_avg !< mean topographic height [Z ~> m] + real, intent(in) :: eta_layer !< height of interface [Z ~> m] real, intent(out) :: w_layer !< frac. open interface width of current layer [nondim] - real, intent(out) :: A_layer !< frac. open face area of current layer [nondim] + real, intent(out) :: A_layer !< frac. open face area of current layer [Z ~> m] + real, intent(in) :: Z_ref !< reference value for geometric height fields [Z ~> m] + real, intent(in) :: Z_to_m !< a unit conversion factor from units of Z to m !local variables - real m, a, & !convenience constant for fit [nondim] - zeta, & !normalized vertical coordinate [nondim] - psi, & !fractional width of layer between Dmin and Dmax [nondim] - psi_int !integral of psi from 0 to zeta + real Dmin, Dmax, Davg, & !copies of input topographic heights stored in [m] + etam, & !copy of eta_layer stored in [m] + m, a, & !convenience constant for fit [nondim] + zeta, & !normalized vertical coordinate [nondim] + psi, & !fractional width of layer between Dmin and Dmax [nondim] + psi_int !integral of psi from 0 to zeta + + Dmin = Z_to_m*(D_min - Z_ref); Dmax = Z_to_m*(D_max - Z_ref) + Davg = Z_to_m*(D_avg - Z_ref); etam = Z_to_m*(eta_layer - Z_ref) !three parameter fit from Adcroft 2013 - m = (D_avg - D_min)/(D_max - D_min) + m = (Davg - Dmin)/(Dmax - Dmin) a = (1. - m)/m - zeta = (eta_layer - D_min)/(D_max - D_min) + zeta = (etam - Dmin)/(Dmax - Dmin) - if (eta_layer <= D_min) then + if (etam <= Dmin) then w_layer = 0.0 A_layer = 0.0 - elseif (eta_layer >= D_max) then + elseif (etam >= Dmax) then w_layer = 1.0 - A_layer = eta_layer - D_avg + A_layer = etam - Davg else if (m < 0.5) then psi = zeta**(1./a) @@ -163,7 +166,7 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) psi_int = zeta - m + m*((1-zeta)**(1/m)) endif w_layer = psi - A_layer = (D_max - D_min)*psi_int + A_layer = (Dmax - Dmin)*psi_int endif diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index efa3b02b2b..5b4f34a87e 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -336,6 +336,14 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) call rotate_array_pair(G_in%areaCu, G_in%areaCv, turns, G%areaCu, G%areaCv) call rotate_array_pair(G_in%IareaCu, G_in%IareaCv, turns, G%IareaCu, G%IareaCv) + call rotate_array_pair(G_in%porous_DminU, G_in%porous_DminV, & + turns, G%porous_DminU, G%porous_DminV) + call rotate_array_pair(G_in%porous_DmaxU, G_in%porous_DmaxV, & + turns, G%porous_DmaxU, G%porous_DmaxV) + call rotate_array_pair(G_in%porous_DavgU, G_in%porous_DavgV, & + turns, G%porous_DavgU, G%porous_DavgV) + + ! Vertex point call rotate_array(G_in%geoLonBu, turns, G%geoLonBu) call rotate_array(G_in%geoLatBu, turns, G%geoLatBu) From 04d932e271741dd05a4045c1acae6da9ece514a9 Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Wed, 24 Nov 2021 17:35:16 -0500 Subject: [PATCH 4/8] Remove Z_to_m rescaling in porous module --- src/core/MOM_porous_barriers.F90 | 17 +++++------------ 1 file changed, 5 insertions(+), 12 deletions(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 34af04a7e4..cdaf566342 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -68,14 +68,14 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) if (G%porous_DavgU(I,j) < 0.) then do K = nk+1,1,-1 eta_s = max(eta(I,j,K), eta(I+1,j,K)) !take shallower layer height - if (US%Z_to_m*eta_s <= (US%Z_to_m*G%porous_DminU(I,j))) then + if (eta_s <= G%porous_DminU(I,j)) then pbv%por_layer_widthU(I,j,K) = 0.0 A_layer_prev = 0.0 if (K < nk+1) then pbv%por_face_areaU(I,j,k) = 0.0; endif else call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), & - G%porous_DavgU(I,j), eta_s, w_layer, A_layer, G%Z_ref, US%Z_to_m) + G%porous_DavgU(I,j), eta_s, w_layer, A_layer) pbv%por_layer_widthU(I,j,K) = w_layer if (k <= nk) then if ((eta_s - eta_prev) > 0.0) then @@ -95,14 +95,14 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) if (G%porous_DavgV(i,J) < 0.) then do K = nk+1,1,-1 eta_s = max(eta(i,J,K), eta(i,J+1,K)) !take shallower layer height - if (US%Z_to_m*eta_s <= US%Z_to_m*G%porous_DminV(i,J)) then + if (eta_s <= G%porous_DminV(i,J)) then pbv%por_layer_widthV(i,J,K) = 0.0 A_layer_prev = 0.0 if (K < nk+1) then pbv%por_face_areaV(i,J,k) = 0.0; endif else call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), & - G%porous_DavgV(i,J), eta_s, w_layer, A_layer, G%Z_ref, US%Z_to_m) + G%porous_DavgV(i,J), eta_s, w_layer, A_layer) pbv%por_layer_widthV(i,J,K) = w_layer if (k <= nk) then if ((eta_s - eta_prev) > 0.0) then @@ -129,19 +129,12 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer, Z_re real, intent(in) :: eta_layer !< height of interface [Z ~> m] real, intent(out) :: w_layer !< frac. open interface width of current layer [nondim] real, intent(out) :: A_layer !< frac. open face area of current layer [Z ~> m] - real, intent(in) :: Z_ref !< reference value for geometric height fields [Z ~> m] - real, intent(in) :: Z_to_m !< a unit conversion factor from units of Z to m !local variables - real Dmin, Dmax, Davg, & !copies of input topographic heights stored in [m] - etam, & !copy of eta_layer stored in [m] - m, a, & !convenience constant for fit [nondim] + real m, a, & !convenience constant for fit [nondim] zeta, & !normalized vertical coordinate [nondim] psi, & !fractional width of layer between Dmin and Dmax [nondim] psi_int !integral of psi from 0 to zeta - Dmin = Z_to_m*(D_min - Z_ref); Dmax = Z_to_m*(D_max - Z_ref) - Davg = Z_to_m*(D_avg - Z_ref); etam = Z_to_m*(eta_layer - Z_ref) - !three parameter fit from Adcroft 2013 m = (Davg - Dmin)/(Dmax - Dmin) a = (1. - m)/m From a98c3ef0d85d02b0e053d494e82fbabfebc7c559 Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Wed, 24 Nov 2021 17:39:55 -0500 Subject: [PATCH 5/8] small bug fix --- src/core/MOM_porous_barriers.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index cdaf566342..f3444e6476 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -121,7 +121,7 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) end subroutine por_widths !> subroutine to calculate the profile fit for a single layer in a column -subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer, Z_ref, Z_to_m) +subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) real, intent(in) :: D_min !< minimum topographic height [Z ~> m] real, intent(in) :: D_max !< maximum topographic height [Z ~> m] From 66e1493affbf24ca2334f6b0de448801501c11f1 Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Wed, 24 Nov 2021 17:41:46 -0500 Subject: [PATCH 6/8] small bug fix 2 --- src/core/MOM_porous_barriers.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index f3444e6476..4f4fd087bb 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -80,7 +80,7 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) if (k <= nk) then if ((eta_s - eta_prev) > 0.0) then pbv%por_face_areaU(I,j,k) = (A_layer - A_layer_prev)/& - (US%Z_to_m*(eta_s-eta_prev)) + (eta_s-eta_prev) else pbv%por_face_areaU(I,j,k) = 0.0; endif endif @@ -107,7 +107,7 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) if (k <= nk) then if ((eta_s - eta_prev) > 0.0) then pbv%por_face_areaV(i,J,k) = (A_layer - A_layer_prev)/& - (US%Z_to_m*(eta_s-eta_prev)) + (eta_s-eta_prev) else pbv%por_face_areaU(I,j,k) = 0.0; endif endif From 62786c40014da8463afb35421dc7315815624184 Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Wed, 24 Nov 2021 18:22:16 -0500 Subject: [PATCH 7/8] fix merge issue --- src/core/MOM_porous_barriers.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index 4f4fd087bb..b20fb993e0 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -132,21 +132,21 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) !local variables real m, a, & !convenience constant for fit [nondim] zeta, & !normalized vertical coordinate [nondim] - psi, & !fractional width of layer between Dmin and Dmax [nondim] + psi, & !fractional width of layer between D_min and D_max [nondim] psi_int !integral of psi from 0 to zeta !three parameter fit from Adcroft 2013 - m = (Davg - Dmin)/(Dmax - Dmin) + m = (D_avg - D_min)/(D_max - D_min) a = (1. - m)/m - zeta = (etam - Dmin)/(Dmax - Dmin) + zeta = (eta_layer - D_min)/(D_max - D_min) - if (etam <= Dmin) then + if (eta_layer <= D_min) then w_layer = 0.0 A_layer = 0.0 - elseif (etam >= Dmax) then + elseif (eta_layer >= D_max) then w_layer = 1.0 - A_layer = etam - Davg + A_layer = eta_layer - D_avg else if (m < 0.5) then psi = zeta**(1./a) @@ -159,7 +159,7 @@ subroutine calc_por_layer(D_min, D_max, D_avg, eta_layer, w_layer, A_layer) psi_int = zeta - m + m*((1-zeta)**(1/m)) endif w_layer = psi - A_layer = (Dmax - Dmin)*psi_int + A_layer = (D_max - D_min)*psi_int endif From ec84be587f53ec7a5241952aefcd4b4d41515288 Mon Sep 17 00:00:00 2001 From: sditkovsky <70655988+sditkovsky@users.noreply.github.com> Date: Mon, 29 Nov 2021 13:59:47 -0500 Subject: [PATCH 8/8] change conversion H_to_MKS to H_to_m for u/veffA * also small fix to porous module. invert order of do loops --- src/core/MOM_dynamics_split_RK2.F90 | 4 ++-- src/core/MOM_dynamics_unsplit.F90 | 4 ++-- src/core/MOM_dynamics_unsplit_RK2.F90 | 4 ++-- src/core/MOM_porous_barriers.F90 | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index e3f30417ef..7cf42319e4 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1542,10 +1542,10 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & - 'Effective U-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + 'Effective U-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & y_cell_method='sum', v_extensive = .true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & - 'Effective V-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + 'Effective V-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & x_cell_method='sum', v_extensive = .true.) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index e52085c2d5..99d092a7f6 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -707,10 +707,10 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & - 'Effective U Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + 'Effective U Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & y_cell_method='sum', v_extensive = .true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & - 'Effective V Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + 'Effective V Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & x_cell_method='sum', v_extensive = .true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index f8112fe4cd..9730fb22fd 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -668,10 +668,10 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & - 'Effective U-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + 'Effective U-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & y_cell_method='sum', v_extensive = .true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & - 'Effective V-Face Area', 'm^2', conversion = GV%H_to_MKS*US%L_to_m, & + 'Effective V-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & x_cell_method='sum', v_extensive = .true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index b20fb993e0..d23509d5f6 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -91,7 +91,7 @@ subroutine por_widths(h, tv, G, GV, US, eta, pbv, eta_bt, halo_size, eta_to_m) endif enddo; enddo - do i=isd,ied; do J=JsdB,JedB + do J=JsdB,JedB; do i=isd,ied if (G%porous_DavgV(i,J) < 0.) then do K = nk+1,1,-1 eta_s = max(eta(i,J,K), eta(i,J+1,K)) !take shallower layer height