diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index aff4860a21..25ed1b1f6e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1013,7 +1013,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & 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, & + call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, US, & CS%set_visc_CSp, symmetrize=.true.) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM)") @@ -2204,7 +2204,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & allocate(CS%tv%frazil(isd:ied,jsd:jed)) ; CS%tv%frazil(:,:) = 0.0 endif if (bound_salinity) then - allocate(CS%tv%salt_deficit(isd:ied,jsd:jed)) ; CS%tv%salt_deficit(:,:)=0.0 + allocate(CS%tv%salt_deficit(isd:ied,jsd:jed)) ; CS%tv%salt_deficit(:,:) = 0.0 endif if (bulkmixedlayer .or. use_temperature) then @@ -2369,20 +2369,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (associated(sponge_in_CSp)) then ! TODO: Implementation and testing of non-ALE spong rotation - call MOM_error(FATAL, "Index rotation of non-ALE sponge is not yet " & - // "implemented.") + call MOM_error(FATAL, "Index rotation of non-ALE sponge is not yet implemented.") endif if (associated(ALE_sponge_in_CSp)) then - call rotate_ALE_sponge(ALE_sponge_in_CSp, G_in, CS%ALE_sponge_CSp, G, & - turns, param_file) - call update_ALE_sponge_field(CS%ALE_sponge_CSp, T_in, CS%T) - call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, CS%S) + call rotate_ALE_sponge(ALE_sponge_in_CSp, G_in, CS%ALE_sponge_CSp, G, turns, param_file) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, T_in, G, GV, CS%T) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, G, GV, CS%S) endif if (associated(OBC_in)) & - call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, & - CS%OBC) + call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, CS%OBC) deallocate(u_in) deallocate(v_in) @@ -2427,7 +2424,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & G => CS%G if (CS%debug .or. CS%G%symmetric) then call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.) - else ; CS%G%Domain_aux => CS%G%Domain ;endif + else ; CS%G%Domain_aux => CS%G%Domain ; endif G%ke = GV%ke endif diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 891d6b3ea7..8d012377a5 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -607,15 +607,16 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ISS%tflux_ocn(i,j) = 0.0 endif -! haline_driving(:,:) = sfc_state%sss(i,j) - Sbdry(i,j) +! haline_driving(i,j) = sfc_state%sss(i,j) - Sbdry(i,j) enddo ! i-loop enddo ! j-loop - ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1] - fluxes%iceshelf_melt(:,:) = ISS%water_flux(:,:) * CS%flux_factor do j=js,je ; do i=is,ie + ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1] + fluxes%iceshelf_melt(i,j) = ISS%water_flux(i,j) * CS%flux_factor + if ((sfc_state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then @@ -653,11 +654,10 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ISS%water_flux(i,j) = 0.0 fluxes%iceshelf_melt(i,j) = 0.0 endif ! area_shelf_h - enddo ; enddo ! i- and j-loops - ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags. - mass_flux(:,:) = 0.0 - mass_flux(:,:) = ISS%water_flux(:,:) * ISS%area_shelf_h(:,:) + ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags. + mass_flux(i,j) = ISS%water_flux(i,j) * ISS%area_shelf_h(i,j) + enddo ; enddo ! i- and j-loops if (CS%active_shelf_dynamics .or. CS%override_shelf_movement) then call cpu_clock_begin(id_clock_pass) @@ -690,7 +690,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well.. ! when we decide on how to do it call update_ice_shelf(CS%dCS, ISS, G, US, US%s_to_T*time_step, Time, & - sfc_state%ocean_mass(:,:), coupled_GL) + sfc_state%ocean_mass, coupled_GL) endif @@ -1735,7 +1735,9 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) call time_interp_external(CS%id_read_mass, Time, ISS%mass_shelf) ! This should only be done if time_interp_external did an update. - ISS%mass_shelf(:,:) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(:,:) ! Rescale after time_interp + do j=js,je ; do i=is,ie + ISS%mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(i,j) ! Rescale after time_interp + enddo ; enddo do j=js,je ; do i=is,ie ISS%area_shelf_h(i,j) = 0.0 diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index dda892dc3e..6145fb1dce 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -558,7 +558,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! Output 2-D energy loss (summed over angles) for each freq and mode do m=1,CS%NMode ; do fr=1,CS%Nfreq if (CS%id_itidal_loss_mode(fr,m) > 0 .or. CS%id_allprocesses_loss_mode(fr,m) > 0) then - itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) + itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together do a=1,CS%nAngle ; do j=js,je ; do i=is,ie itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + CS%TKE_itidal_loss(i,j,a,fr,m) @@ -886,12 +886,16 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) !! across angles [R Z3 T-2 ~> J m-2]. ! Local variables real :: flux - real :: u_ang - real :: Angle_size - real :: I_Angle_size - real :: I_dt + real :: u_ang ! Angular propagation speed [Rad T-1 ~> Rad s-1] + real :: Angle_size ! The size of each orientation wedge in radians [Rad] + real :: I_Angle_size ! The inverse of the the orientation wedges [Rad-1] + real :: I_dt ! The inverse of the timestep [T-1 ~> s-1] + real :: aR, aL ! Left and right edge estimates of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real :: dMx, dMn + real :: Ep, Ec, Em ! Mean angular energy density for three successive wedges in angular + ! orientation [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real :: dA, mA, a6 ! Difference, mean, and curvature of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] integer :: a - real :: aR, aL, dMx, dMn, Ep, Ec, Em, dA, mA, a6 I_dt = 1 / dt Angle_size = (8.0*atan(1.0)) / (real(NAngle)) @@ -902,9 +906,12 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) u_ang = CFL_ang(A)*Angle_size*I_dt if (u_ang >= 0.0) then ! Implementation of PPM-H3 - Ep = En2d(a+1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Ec = En2d(a) *I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Em = En2d(a-1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) + ! Convert wedge-integrated energy density into angular energy densities for three successive + ! wedges around the source wedge for this flux [R Z3 T-2 rad-1 ~> J m-2 rad-1]. + Ep = En2d(a+1)*I_Angle_size + Ec = En2d(a) *I_Angle_size + Em = En2d(a-1)*I_Angle_size + ! Calculate and bound edge values of energy density. aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate @@ -918,17 +925,21 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) aR = 3.*Ec - 2.*aL !? endif a6 = 6.*Ec - 3. * (aR + aL) ! Curvature - ! CALCULATE FLUX RATE (Jm-2/s) + ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - !flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! CALCULATE AMOUNT FLUXED (Jm-2) + ! The following expression copied from tracer_advect is equivalent. + ! flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) + ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux else ! Implementation of PPM-H3 - Ep = En2d(a+2)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Ec = En2d(a+1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Em = En2d(a) *I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) + ! Convert wedge-integrated energy density into angular energy densities for three successive + ! wedges around the source wedge for this flux [R Z3 T-2 rad-1 ~> J m-2 rad-1]. + Ep = En2d(a+2)*I_Angle_size + Ec = En2d(a+1)*I_Angle_size + Em = En2d(a) *I_Angle_size + ! Calculate and bound edge values of energy density. aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate @@ -942,10 +953,12 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) aR = 3.*Ec - 2.*aL endif a6 = 6.*Ec - 3. * (aR + aL) ! Curvature - ! CALCULATE FLUX RATE (Jm-2/s) + ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] + !### This expression is wrong, because it was just copied from above. The correct one is below flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - !flux = u_ang*( aL + 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! CALCULATE AMOUNT FLUXED (Jm-2) + ! This is the correct expression; note that CFL_ang is negative here, so it looks a bit odd. + !flux = u_ang*( aL - 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. + 2./3. * CFL_ang(A) ) ) ) + ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux endif @@ -1014,7 +1027,7 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle) ! FIND AVERAGE GROUP VELOCITY (SPEED) AT CELL CORNERS ! NOTE: THIS HAS NOT BE ADAPTED FOR REFLECTION YET (BDM)!! ! Fix indexing here later - speed(:,:) = 0 + speed(:,:) = 0.0 do J=jsh-1,jeh ; do I=ish-1,ieh f2 = G%CoriolisBu(I,J)**2 speed(I,J) = 0.25*(cn(i,j) + cn(i+1,j) + cn(i+1,j+1) + cn(i,j+1)) * & @@ -1058,21 +1071,21 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle) ! Apply propagation in x-direction (reflection included) LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh - call propagate_x(En(:,:,:), speed_x, Cgx_av(:), dCgx(:), dt, G, US, CS%nAngle, CS, LB) + call propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, CS%nAngle, CS, LB) ! Check for energy conservation on computational domain (for debugging) - !call sum_En(G,CS,En(:,:,:),'post-propagate_x') + !call sum_En(G, CS, En, 'post-propagate_x') ! Update halos - call pass_var(En(:,:,:),G%domain) + call pass_var(En, G%domain) ! Apply propagation in y-direction (reflection included) ! LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie ! Use if no teleport LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh - call propagate_y(En(:,:,:), speed_y, Cgy_av(:), dCgy(:), dt, G, US, CS%nAngle, CS, LB) + call propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, CS%nAngle, CS, LB) ! Check for energy conservation on computational domain (for debugging) - !call sum_En(G,CS,En(:,:,:),'post-propagate_y') + !call sum_En(G, CS, En, 'post-propagate_y') endif end subroutine propagate @@ -1084,7 +1097,7 @@ subroutine propagate_corner_spread(En, energized_wedge, NAngle, speed, dt, G, CS type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, dimension(G%isd:G%ied,G%jsd:G%jed), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed), & intent(in) :: speed !< The magnitude of the group velocity at the cell !! corner points [L T-1 ~> m s-1]. @@ -1351,7 +1364,7 @@ subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB) !! discretized wave energy spectrum. real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), & intent(in) :: speed_x !< The magnitude of the group velocity at the !! Cu points [L T-1 ~> m s-1]. @@ -1404,18 +1417,18 @@ subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB) enddo ! a-loop ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected - ! and will eventually propagate out of cell. (Thid code only reflects if En > 0) - call reflect(Fdt_m(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_m(:,:,:), Nangle, CS, G, LB) - call reflect(Fdt_p(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_p(:,:,:), Nangle, CS, G, LB) + ! and will eventually propagate out of cell. (This code only reflects if En > 0.) + call reflect(Fdt_m, Nangle, CS, G, LB) + call teleport(Fdt_m, Nangle, CS, G, LB) + call reflect(Fdt_p, Nangle, CS, G, LB) + call teleport(Fdt_p, Nangle, CS, G, LB) - ! Update reflected energy (Jm-2) - do j=jsh,jeh ; do i=ish,ieh + ! Update reflected energy [R Z3 T-2 ~> J m-2] + do a=1,Nangle ; do j=jsh,jeh ; do i=ish,ieh ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging ! call MOM_error(FATAL, "propagate_x: OutFlux>Available") - En(i,j,:) = En(i,j,:) + G%IareaT(i,j)*(Fdt_m(i,j,:) + Fdt_p(i,j,:)) - enddo ; enddo + En(i,j,a) = En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a)) + enddo ; enddo ; enddo end subroutine propagate_x @@ -1426,7 +1439,7 @@ subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB) !! discretized wave energy spectrum. real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%isd:G%ied,G%JsdB:G%JedB), & intent(in) :: speed_y !< The magnitude of the group velocity at the !! Cv points [L T-1 ~> m s-1]. @@ -1486,13 +1499,13 @@ subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB) enddo ! a-loop ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected - ! and will eventually propagate out of cell. (Thid code only reflects if En > 0) - call reflect(Fdt_m(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_m(:,:,:), Nangle, CS, G, LB) - call reflect(Fdt_p(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_p(:,:,:), Nangle, CS, G, LB) + ! and will eventually propagate out of cell. (This code only reflects if En > 0.) + call reflect(Fdt_m, Nangle, CS, G, LB) + call teleport(Fdt_m, Nangle, CS, G, LB) + call reflect(Fdt_p, Nangle, CS, G, LB) + call teleport(Fdt_p, Nangle, CS, G, LB) - ! Update reflected energy (Jm-2) + ! Update reflected energy [R Z3 T-2 ~> J m-2] do a=1,Nangle ; do j=jsh,jeh ; do i=ish,ieh ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging ! call MOM_error(FATAL, "propagate_y: OutFlux>Available", .true.) @@ -1521,8 +1534,7 @@ subroutine zonal_flux_En(u, h, hL, hR, uh, dt, G, US, j, ish, ieh, vol_CFL) !! the cell areas when estimating the CFL number. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]. - real :: curv_3 ! A measure of the thickness curvature over a grid length, - ! with the same units as h_in. + real :: curv_3 ! A measure of the energy density curvature over a grid length [R Z3 T-2 ~> J m-2] integer :: i do I=ish-1,ieh @@ -1566,8 +1578,7 @@ subroutine merid_flux_En(v, h, hL, hR, vh, dt, G, US, J, ish, ieh, vol_CFL) !! the CFL number. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]. - real :: curv_3 ! A measure of the thickness curvature over a grid length, - ! with the same units as h_in. + real :: curv_3 ! A measure of the energy density curvature over a grid length [R Z3 T-2 ~> J m-2] integer :: i do i=ish,ieh @@ -1603,18 +1614,18 @@ subroutine reflect(En, NAngle, CS, G, LB) type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds. ! Local variables real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c - ! angle of boudary wrt equator + ! angle of boundary wrt equator [rad] real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl ! fraction of wave energy reflected - ! values should collocate with angle_c + ! values should collocate with angle_c [nondim] logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge ! tags of cells with double reflection - real :: TwoPi ! 2*pi - real :: Angle_size ! size of beam wedge (rad) - real :: angle_wall ! angle of coast/ridge/shelf wrt equator - real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator - real :: angle_r ! angle of reflected ray wrt equator + real :: TwoPi ! 2*pi = 6.2831853... [nondim] + real :: Angle_size ! size of beam wedge [rad] + real :: angle_wall ! angle of coast/ridge/shelf wrt equator [rad] + real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator [rad] + real :: angle_r ! angle of reflected ray wrt equator [rad] real, dimension(1:Nangle) :: En_reflected integer :: i, j, a, a_r, na !integer :: isd, ied, jsd, jed ! start and end local indices on data domain @@ -1623,7 +1634,6 @@ subroutine reflect(En, NAngle, CS, G, LB) ! (values exclude halos) integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain ! leaving out outdated halo points (march in) - integer :: id_g, jd_g ! global (decomp-invar) indices !isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec @@ -1643,59 +1653,54 @@ subroutine reflect(En, NAngle, CS, G, LB) ridge = CS%refl_dbl En_reflected(:) = 0.0 - !do j=jsc-1,jec+1 - do j=jsh,jeh - !do i=isc-1,iec+1 - do i=ish,ieh - ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset - ! redistribute energy in angular space if ray will hit boundary - ! i.e., if energy is in a reflecting cell - if (angle_c(i,j) /= CS%nullangle) then - do a=1,NAngle - if (En(i,j,a) > 0.0) then - ! if ray is incident, keep specified boundary angle - if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then - angle_wall = angle_c(i,j) - ! if ray is not incident but in ridge cell, use complementary angle - elseif (ridge(i,j)) then - angle_wall = angle_c(i,j) + 0.5*TwoPi - if (angle_wall > TwoPi) then - angle_wall = angle_wall - TwoPi*floor(abs(angle_wall)/TwoPi) - elseif (angle_wall < 0.0) then - angle_wall = angle_wall + TwoPi*ceiling(abs(angle_wall)/TwoPi) - endif - ! if ray is not incident and not in a ridge cell, keep specified angle - else - angle_wall = angle_c(i,j) - endif - ! do reflection - if (sin(angle_i(a) - angle_wall) >= 0.0) then - angle_r = 2.0*angle_wall - angle_i(a) - if (angle_r > TwoPi) then - angle_r = angle_r - TwoPi*floor(abs(angle_r)/TwoPi) - elseif (angle_r < 0.0) then - angle_r = angle_r + TwoPi*ceiling(abs(angle_r)/TwoPi) - endif - a_r = nint(angle_r/Angle_size) + 1 - do while (a_r > Nangle) ; a_r = a_r - Nangle ; enddo - if (a /= a_r) then - En_reflected(a_r) = part_refl(i,j)*En(i,j,a) - En(i,j,a) = (1.0-part_refl(i,j))*En(i,j,a) - endif - endif + do j=jsh,jeh ; do i=ish,ieh + ! redistribute energy in angular space if ray will hit boundary + ! i.e., if energy is in a reflecting cell + if (angle_c(i,j) /= CS%nullangle) then + do a=1,NAngle ; if (En(i,j,a) > 0.0) then + if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then + ! if ray is incident, keep specified boundary angle + angle_wall = angle_c(i,j) + elseif (ridge(i,j)) then + ! if ray is not incident but in ridge cell, use complementary angle + angle_wall = angle_c(i,j) + 0.5*TwoPi + if (angle_wall > TwoPi) then + angle_wall = angle_wall - TwoPi*floor(abs(angle_wall)/TwoPi) + elseif (angle_wall < 0.0) then + angle_wall = angle_wall + TwoPi*ceiling(abs(angle_wall)/TwoPi) endif - enddo ! a-loop - En(i,j,:) = En(i,j,:) + En_reflected(:) - En_reflected(:) = 0.0 - endif - enddo ! i-loop - enddo ! j-loop + else + ! if ray is not incident and not in a ridge cell, keep specified angle + angle_wall = angle_c(i,j) + endif + + ! do reflection + if (sin(angle_i(a) - angle_wall) >= 0.0) then + angle_r = 2.0*angle_wall - angle_i(a) + if (angle_r > TwoPi) then + angle_r = angle_r - TwoPi*floor(abs(angle_r)/TwoPi) + elseif (angle_r < 0.0) then + angle_r = angle_r + TwoPi*ceiling(abs(angle_r)/TwoPi) + endif + a_r = nint(angle_r/Angle_size) + 1 + do while (a_r > Nangle) ; a_r = a_r - Nangle ; enddo + if (a /= a_r) then + En_reflected(a_r) = part_refl(i,j)*En(i,j,a) + En(i,j,a) = (1.0-part_refl(i,j))*En(i,j,a) + endif + endif + endif ; enddo ! a-loop + do a=1,NAngle + En(i,j,a) = En(i,j,a) + En_reflected(a) + En_reflected(a) = 0.0 + enddo ! a-loop + endif + enddo ; enddo ! i- and j-loops ! Check to make sure no energy gets onto land (only run for debugging) ! do a=1,NAngle ; do j=jsc,jec ; do i=isc,iec ! if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then - ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset - ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',id_g, 'jg_g=',jd_g + ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',i+G%idg_offset, 'jg_g=',j+G%jdg_offset ! call MOM_error(FATAL, "reflect: Energy detected out of bounds: "//trim(mesg), .true.) ! endif ! enddo ; enddo ; enddo @@ -1717,17 +1722,17 @@ subroutine teleport(En, NAngle, CS, G, LB) type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds. ! Local variables real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c - ! angle of boudary wrt equator + ! angle of boundary wrt equator [rad] real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl ! fraction of wave energy reflected - ! values should collocate with angle_c + ! values should collocate with angle_c [nondim] logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell ! flag for partial reflection logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge - ! tags of cells with double reflection - real :: TwoPi ! size of beam wedge (rad) - real :: Angle_size ! size of beam wedge (rad) - real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator + ! tags of cells with double reflection + real :: TwoPi ! 2*pi = 6.2831853... [nondim] + real :: Angle_size ! size of beam wedge [rad] + real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator [rad] real, dimension(1:NAngle) :: cos_angle, sin_angle real :: En_tele ! energy to be "teleported" [R Z3 T-2 ~> J m-2] character(len=160) :: mesg ! The text of an error message @@ -2295,8 +2300,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) CS%TKE_itidal_loss(:,:,:,:,:) = 0.0 allocate(CS%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode)) CS%TKE_Froude_loss(:,:,:,:,:) = 0.0 - allocate(CS%tot_leak_loss(isd:ied,jsd:jed)) ; CS%tot_leak_loss(:,:) = 0.0 - allocate(CS%tot_quad_loss(isd:ied,jsd:jed) ) ; CS%tot_quad_loss(:,:) = 0.0 + allocate(CS%tot_leak_loss(isd:ied,jsd:jed)) ; CS%tot_leak_loss(:,:) = 0.0 + allocate(CS%tot_quad_loss(isd:ied,jsd:jed) ) ; CS%tot_quad_loss(:,:) = 0.0 allocate(CS%tot_itidal_loss(isd:ied,jsd:jed)) ; CS%tot_itidal_loss(:,:) = 0.0 allocate(CS%tot_Froude_loss(isd:ied,jsd:jed)) ; CS%tot_Froude_loss(:,:) = 0.0 diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index b791535ed1..fe1ccab53d 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -171,7 +171,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v character(len=10) :: remapScheme if (associated(CS)) then - call MOM_error(WARNING, "initialize_sponge called with an associated "// & + call MOM_error(WARNING, "initialize_ALE_sponge_fixed called with an associated "// & "control structure.") return endif @@ -260,14 +260,14 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if (CS%sponge_uv) then - allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)); data_hu(:,:,:)=0.0 - allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)); data_hv(:,:,:)=0.0 - allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)); Iresttime_u(:,:)=0.0 - allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)); Iresttime_v(:,:)=0.0 + allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)) ; data_hu(:,:,:) = 0.0 + allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)) ; data_hv(:,:,:) = 0.0 + allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 + allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB + do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB data_hu(I,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:)) Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & @@ -276,9 +276,9 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if (CS%num_col_u > 0) then - allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u = 0.0 - allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u = 0 - allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u = 0 + allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u(:) = 0.0 + allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u(:) = 0 + allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u(:) = 0 ! pass indices, restoring time to the CS structure col = 1 @@ -286,7 +286,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = i ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo @@ -323,7 +323,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo @@ -415,7 +415,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) character(len=10) :: remapScheme if (associated(CS)) then - call MOM_error(WARNING, "initialize_sponge called with an associated "// & + call MOM_error(WARNING, "initialize_ALE_sponge_varying called with an associated "// & "control structure.") return endif @@ -486,8 +486,8 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, & "The total number of columns where sponges are applied at h points.") if (CS%sponge_uv) then - allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)); Iresttime_u(:,:)=0.0 - allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)); Iresttime_v(:,:)=0.0 + allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 + allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB @@ -578,7 +578,7 @@ subroutine set_up_ALE_sponge_field_fixed(sp_val, G, f_ptr, CS) if (CS%fldno > MAX_FIELDS_) then write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & &the number of fields to be damped in the call to & - &initialize_sponge." )') CS%fldno + &initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif @@ -605,8 +605,8 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(in) :: G !< Grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & target, intent(in) :: f_ptr !< Pointer to the field to be damped (in). type(ALE_sponge_CS), pointer :: CS !< Sponge control structure (in/out). @@ -634,7 +634,7 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, if (CS%fldno > MAX_FIELDS_) then write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & &the number of fields to be damped in the call to & - &initialize_sponge." )') CS%fldno + &initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif ! get a unique time interp id for this field. If sponge data is ongrid, then setup @@ -788,11 +788,11 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]. type(ALE_sponge_CS), pointer :: CS !< A pointer to the control structure for this module - !! that is set by a previous call to initialize_sponge (in). + !! that is set by a previous call to initialize_ALE_sponge (in). type(time_type), optional, intent(in) :: Time !< The current model date real :: damp ! The timestep times the local damping coefficient [nondim]. @@ -833,8 +833,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) nz_data = CS%Ref_val(m)%nz_data allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) - sp_val(:,:,:)=0.0 - mask_z(:,:,:)=0.0 + sp_val(:,:,:) = 0.0 + mask_z(:,:,:) = 0.0 call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & z_edges_in, missing_value, .true., .false., .false., & spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & @@ -1003,17 +1003,20 @@ end subroutine apply_ALE_sponge !> Rotate the ALE sponge fields from the input to the model index map. subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) - type(ALE_sponge_CS), intent(in) :: sponge_in - type(ocean_grid_type), intent(in) :: G_in - type(ALE_sponge_CS), pointer :: sponge - type(ocean_grid_type), intent(in) :: G - integer, intent(in) :: turns - type(param_file_type), intent(in) :: param_file + type(ALE_sponge_CS), intent(in) :: sponge_in !< The control structure for this module with the + !! original grid rotation + type(ocean_grid_type), intent(in) :: G_in !< The ocean's grid structure with the original rotation. + type(ALE_sponge_CS), pointer :: sponge !< A pointer to the control that will be set up with + !! the new grid rotation + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure with the new rotation. + integer, intent(in) :: turns !< The number of 90-degree turns between grids + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. ! First part: Index construction ! 1. Reconstruct Iresttime(:,:) from sponge_in ! 2. rotate Iresttime(:,:) - ! 3. Call initialize_sponge using new grid and rotated Iresttime(:,:) + ! 3. Call initialize_ALE_sponge using new grid and rotated Iresttime(:,:) ! All the index adjustment should follow from the Iresttime rotation real, dimension(:,:), allocatable :: Iresttime_in, Iresttime @@ -1040,15 +1043,13 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) endif ! Re-populate the 2D Iresttime and data_h arrays on the original grid - do c = 1, sponge_in%num_col + do c=1,sponge_in%num_col c_i = sponge_in%col_i(c) c_j = sponge_in%col_j(c) Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) - if (fixed_sponge) then - do k = 1, nz_data - data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) - enddo - endif + if (fixed_sponge) then ; do k=1,nz_data + data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) + enddo ; endif enddo call rotate_array(Iresttime_in, turns, Iresttime) @@ -1076,19 +1077,13 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) allocate(sp_val(G%isd:G%ied, G%jsd:G%jed, nz_data)) endif - do n = 1, sponge_in%fldno + do n=1,sponge_in%fldno ! Assume that tracers are pointers and are remapped in other functions(?) sp_ptr => sponge_in%var(n)%p sp_val_in(:,:,:) = 0.0 - do c = 1, sponge_in%num_col - c_i = sponge_in%col_i(c) - c_j = sponge_in%col_j(c) - if (fixed_sponge) then - do k = 1, nz_data - sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) - enddo - endif - enddo + if (fixed_sponge) then ; do c=1,sponge_in%num_col ; do k=1,nz_data + sp_val_in(sponge_in%col_i(c), sponge_in%col_j(c), k) = sponge_in%Ref_val(n)%p(k,c) + enddo ; enddo ; endif call rotate_array(sp_val_in, turns, sp_val) if (fixed_sponge) then @@ -1144,17 +1139,22 @@ end subroutine rotate_ALE_sponge ! TODO: This function solely exists to replace field pointers in the sponge ! after rotation. This function is part of a temporary solution until ! something more robust is developed. -subroutine update_ALE_sponge_field(sponge, p_old, p_new) - type(ALE_sponge_CS), pointer :: sponge - real, dimension(:,:,:), target, intent(in) :: p_old - real, dimension(:,:,:), target, intent(in) :: p_new +subroutine update_ALE_sponge_field(sponge, p_old, G, GV, p_new) + type(ALE_sponge_CS), pointer :: sponge !< A pointer to the control structure for this module + !! that is set by a previous call to initialize_ALE_sponge. + real, dimension(:,:,:), & + target, intent(in) :: p_old !< The previous array of target values + type(ocean_grid_type), intent(in) :: G !< The updated ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + target, intent(in) :: p_new !< The new array of target values integer :: n - do n = 1, sponge%fldno - if (associated(sponge%var(n)%p, p_old)) & - sponge%var(n)%p => p_new + do n=1,sponge%fldno + if (associated(sponge%var(n)%p, p_old)) sponge%var(n)%p => p_new enddo + end subroutine update_ALE_sponge_field @@ -1163,7 +1163,7 @@ end subroutine update_ALE_sponge_field !> This subroutine deallocates any memory associated with the ALE_sponge module. subroutine ALE_sponge_end(CS) type(ALE_sponge_CS), pointer :: CS !< A pointer to the control structure that is - !! set by a previous call to initialize_sponge. + !! set by a previous call to initialize_ALE_sponge. integer :: m diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 27b316e144..e7d5bcc476 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -509,8 +509,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, CS) filename) call safe_alloc_ptr(CS%TKE_Niku,is,ie,js,je) ; CS%TKE_Niku(:,:) = 0.0 call MOM_read_data(filename, 'TKE_input', CS%TKE_Niku, G%domain, timelevel=1, & ! ??? timelevel -aja - scale=US%W_m2_to_RZ3_T3) - CS%TKE_Niku(:,:) = Niku_scale * CS%TKE_Niku(:,:) + scale=Niku_scale*US%W_m2_to_RZ3_T3) call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, & "The fraction of the lee wave energy that is dissipated "//& @@ -781,7 +780,7 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k = 1,G%ke+1 + do k=1,G%ke+1 N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) enddo @@ -876,7 +875,7 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) CVmix_tidal_params_user = CS%CVMix_tidal_params) ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k = 1,G%ke+1 + do k=1,G%ke+1 N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) enddo diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index c9bf98f7ff..5503287c50 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -287,7 +287,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G if (.not.associated(CS)) return - melt(:,:) = fluxes%iceshelf_melt + melt(:,:) = fluxes%iceshelf_melt(:,:) ! max. melt mmax = MAXVAL(melt(is:ie,js:je)) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 8b9be533d5..f244931376 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -171,16 +171,14 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) ! for diagnostics if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then - tendency(:,:,:) = 0.0 + tendency(:,:,:) = 0.0 endif - do j = G%jsc-1, G%jec+1 - ! Interpolate state to interface - do i = G%isc-1, G%iec+1 - call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & - ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) - enddo - enddo + ! Interpolate state to interface + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & + ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) + enddo ; enddo ! Diffusive fluxes in the i-direction uFlx(:,:,:) = 0. vFlx(:,:,:) = 0. @@ -253,41 +251,41 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vFlx*Idt, CS%diag) if (tracer%id_lbd_dfx_2d>0) then uwork_2d(:,:) = 0. - do k=1,GV%ke; do j=G%jsc,G%jec; do I=G%isc-1,G%iec + do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec uwork_2d(I,j) = uwork_2d(I,j) + (uFlx(I,j,k) * Idt) - enddo; enddo; enddo + enddo ; enddo ; enddo call post_data(tracer%id_lbd_dfx_2d, uwork_2d, CS%diag) endif if (tracer%id_lbd_dfy_2d>0) then vwork_2d(:,:) = 0. - do k=1,GV%ke; do J=G%jsc-1,G%jec; do i=G%isc,G%iec + do k=1,GV%ke ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec vwork_2d(i,J) = vwork_2d(i,J) + (vFlx(i,J,k) * Idt) - enddo; enddo; enddo + enddo ; enddo ; enddo call post_data(tracer%id_lbd_dfy_2d, vwork_2d, CS%diag) endif ! post tendency of tracer content if (tracer%id_lbdxy_cont > 0) then - call post_data(tracer%id_lbdxy_cont, tendency(:,:,:), CS%diag) + call post_data(tracer%id_lbdxy_cont, tendency, CS%diag) endif ! post depth summed tendency for tracer content if (tracer%id_lbdxy_cont_2d > 0) then tendency_2d(:,:) = 0. - do j = G%jsc,G%jec ; do i = G%isc,G%iec - do k = 1, GV%ke + do j=G%jsc,G%jec ; do i=G%isc,G%iec + do k=1,GV%ke tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k) enddo enddo ; enddo - call post_data(tracer%id_lbdxy_cont_2d, tendency_2d(:,:), CS%diag) + call post_data(tracer%id_lbdxy_cont_2d, tendency_2d, CS%diag) endif ! post tendency of tracer concentration; this step must be ! done after posting tracer content tendency, since we alter ! the tendency array and its units. if (tracer%id_lbdxy_conc > 0) then - do k = 1, GV%ke ; do j = G%jsc,G%jec ; do i = G%isc,G%iec + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff ) enddo ; enddo ; enddo call post_data(tracer%id_lbdxy_conc, tendency, CS%diag) @@ -306,9 +304,9 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2] real :: hBLT !< Depth of the boundary layer [H ~> m or kg m-2] real, dimension(nk) :: phi !< Scalar quantity - real, dimension(nk,2) :: ppoly0_E(:,:) !< Edge value of polynomial - real, dimension(nk,deg+1) :: ppoly0_coefs(:,:) !< Coefficients of polynomial - integer :: method !< Remapping scheme to use + real, dimension(nk,2) :: ppoly0_E !< Edge value of polynomial + real, dimension(nk,deg+1) :: ppoly0_coefs !< Coefficients of polynomial + integer :: method !< Remapping scheme to use integer :: k_top !< Index of the first layer within the boundary real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 65f83ecfea..a0f7b23346 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -722,9 +722,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, CS, h_pre, eatr, e ! Add diurnal cycle for shortwave radiation (only used if run in ocean-only mode) if (CS%diurnal_SW .and. CS%read_sw) then - sw(:,:) = fluxes%sw - sw_vis(:,:) = fluxes%sw_vis_dir - sw_nir(:,:) = fluxes%sw_nir_dir + sw(:,:) = fluxes%sw(:,:) + sw_vis(:,:) = fluxes%sw_vis_dir(:,:) + sw_nir(:,:) = fluxes%sw_nir_dir(:,:) call offline_add_diurnal_SW(fluxes, CS%G, Time_start, Time_end) endif @@ -738,9 +738,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, CS, h_pre, eatr, e CS%G, CS%GV, CS%US, CS%tv, CS%optics, CS%tracer_flow_CSp, CS%debug) if (CS%diurnal_SW .and. CS%read_sw) then - fluxes%sw(:,:) = sw - fluxes%sw_vis_dir(:,:) = sw_vis - fluxes%sw_nir_dir(:,:) = sw_nir + fluxes%sw(:,:) = sw(:,:) + fluxes%sw_vis_dir(:,:) = sw_vis(:,:) + fluxes%sw_nir_dir(:,:) = sw_nir(:,:) endif if (CS%debug) then