diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 30b8b271d5..fa350742b6 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -106,6 +106,7 @@ module MOM_ALE public ALE_init public ALE_end public ALE_main +public ALE_main_offline public ALE_offline_tracer_final public ALE_build_grid public ALE_remap_scalar @@ -425,6 +426,61 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt) end subroutine ALE_main +!> Takes care of (1) building a new grid and (2) remapping all variables between +!! the old grid and the new grid. The creation of the new grid can be based +!! on z coordinates, target interface densities, sigma coordinates or any +!! arbitrary coordinate system. +subroutine ALE_main_offline( G, GV, h, tv, Reg, CS, dt) + type(ocean_grid_type), intent(in) :: G !< Ocean grid informations + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after last time step (m or Pa) + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure + type(tracer_registry_type), pointer :: Reg !< Tracer registry structure + type(ALE_CS), pointer :: CS !< Regridding parameters and options + real, optional, intent(in) :: dt !< Time step between calls to ALE_main() + ! Local variables + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step (m or Pa) + integer :: nk, i, j, k, isc, iec, jsc, jec + + nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + + if (CS%show_call_tree) call callTree_enter("ALE_main_offline(), MOM_ALE.F90") + + if (present(dt)) then + call ALE_update_regrid_weights( dt, CS ) + endif + dzRegrid(:,:,:) = 0.0 + + ! Build new grid. The new grid is stored in h_new. The old grid is h. + ! Both are needed for the subsequent remapping of variables. + call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid ) + + call check_grid( G, GV, h, 0. ) + + if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_main)") + + ! Remap all variables from old grid h onto new grid h_new + + call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, -dzRegrid, Reg, & + debug=CS%show_call_tree, dt=dt ) + + if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_main)") + + ! Override old grid with new one. The new grid 'h_new' is built in + ! one of the 'build_...' routines above. +!$OMP parallel do default(none) shared(isc,iec,jsc,jec,nk,h,h_new,CS) + do k = 1,nk + do j = jsc-1,jec+1 ; do i = isc-1,iec+1 + h(i,j,k) = h_new(i,j,k) + enddo ; enddo + enddo + + if (CS%show_call_tree) call callTree_leave("ALE_main()") + if (CS%id_dzRegrid>0 .and. present(dt)) call post_data(CS%id_dzRegrid, dzRegrid, CS%diag) + +end subroutine ALE_main_offline + !> Remaps all tracers from h onto h_target. This is intended to be called when tracers !! are done offline. In the case where transports don't quite conserve, we still want to !! make sure that layer thicknesses offline do not drift too far away from the online model diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f95fbb0ae1..4fbcb5931b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -63,7 +63,7 @@ module MOM use MOM_ALE, only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity use MOM_ALE, only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile use MOM_ALE, only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags -use MOM_ALE, only : ALE_offline_tracer_final +use MOM_ALE, only : ALE_main_offline, ALE_offline_tracer_final use MOM_continuity, only : continuity, continuity_init, continuity_CS use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_init, CoriolisAdv_CS use MOM_diabatic_driver, only : diabatic, diabatic_driver_init, diabatic_CS @@ -131,7 +131,6 @@ module MOM use MOM_offline_transport, only : transport_by_files, next_modulo_time, post_advection_fields use MOM_offline_transport, only : offline_transport_init, register_diags_offline_transport use MOM_offline_transport, only : limit_mass_flux_3d, update_h_horizontal_flux, update_h_vertical_flux -use MOM_offline_transport, only : limit_mass_flux_ordered_3d use MOM_tracer_diabatic, only : applyTracerBoundaryFluxesInOut use time_manager_mod, only : print_date @@ -1502,14 +1501,10 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) ! Read in all fields that might be used this timestep call transport_by_files(G, GV, CS%offline_CSp, h_end, eatr, ebtr, uhtr, vhtr, & khdt_x, khdt_y, temp_old, salt_old, fluxes, CS%use_ALE_algorithm) - if (CS%offline_CSp%id_uhtr_preadv>0) call post_data(CS%offline_CSp%id_uhtr_preadv, uhtr, CS%diag) - if (CS%offline_CSp%id_vhtr_preadv>0) call post_data(CS%offline_CSp%id_vhtr_preadv, vhtr, CS%diag) - if (CS%id_h>0) call post_data(CS%id_h, h_end, CS%diag) do k=1,nz ; do j=jsd,jed ; do i=isd,ied h_pre(i,j,k) = CS%h(i,j,k) enddo ; enddo; enddo - call pass_var(h_pre,G%Domain) @@ -1654,6 +1649,25 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) elseif (CS%use_ALE_algorithm) then + ! Tracers are transported using the stored mass fluxes in the following way + ! 1) Using the layer thicknesses and tracer concentrations from the previous timestep, + ! half of the accumulated vertical mixing (eatr and ebtr) is applied in the call to tracer_column_fns. + ! For tracers whose source/sink terms need dt, this value is set to 1/2 dt_offline + ! 2) Next half of the accumulated surface freshwater fluxes are applied + !! Steps 3, 4, and 5 are iterated + ! 3) Accumulated mass fluxes are used to do horizontal transport. The number of iterations used in + ! advect_tracer is limited to 2 (e.g x->y->x->y). The remaining mass fluxes are stored for later use + ! and resulting layer thicknesses fed into the next step + ! 4) Tracers and the h-grid are regridded and remapped in a call to ALE. This allows for layers which might + ! 'vanish' because of horizontal mass transport to be 'reinflated' + ! 5) Check that transport is done if the remaining mass fluxes equals 0 or if the max number of iterations + ! has been reached + !! END ITERATION + ! 6) Repeat steps 1 and 2 + ! 7) Force a remapping to the stored layer thicknesses that correspond to the snapshot of the online model + ! 8) Reset T/S and h to their stored snapshotted values to prevent model drift + + ! Convert flux rates into explicit mass/height of freshwater flux. Also note, that ! fluxes are halved because diabatic processes are split before and after advection do j=jsd,jed ; do i=isd,ied @@ -1662,7 +1676,8 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) enddo ; enddo zero_3dh(:,:,:)=0.0 - + + ! Copy over the horizontal mass fluxes from the remaining total mass fluxes do k=1,nz ; do j=jsd,jed ; do i=isdB,iedB uhtr_sub(i,j,k) = uhtr(i,j,k) enddo ; enddo ; enddo @@ -1670,6 +1685,14 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) vhtr_sub(i,j,k) = vhtr(i,j,k) enddo ; enddo ; enddo + if(CS%debug) then + call uchksum(uhtr_sub,"uhtr_sub before transport",G%HI) + call vchksum(vhtr_sub,"vhtr_sub before transport",G%HI) + call hchksum(h_pre,"h_pre before transport",G%HI) + endif + + ! Note that here, h_new really shouldn't be used, should double check that any individual tracer + ! does not use h_new call call_tracer_column_fns(h_pre, h_new, eatr*0.5, ebtr*0.5, & fluxes, CS%offline_CSp%dt_offline*0.5, G, GV, CS%tv, & CS%diabatic_CSp%optics, CS%tracer_flow_CSp, CS%debug, & @@ -1677,7 +1700,10 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) minimum_forcing_depth=0.001) call applyTracerBoundaryFluxesInOut(G, GV, zero_3dh, 0.5*dt_offline, fluxes, h_pre, & 0.8, 0.001) - + + if(CS%debug) then + call hchksum(h_pre,"h_pre after 1st diabatic",G%HI) + endif do iter=1,CS%offline_CSp%num_off_iter @@ -1687,18 +1713,27 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) call advect_tracer(h_pre, uhtr_sub, vhtr_sub, CS%OBC, dt_iter, G, GV, & CS%tracer_adv_CSp, CS%tracer_Reg, h_vol, max_iter_in=2, & - uhr_out=uhtr, vhr_out=vhtr, h_out=h_new) + uhr_out=uhtr, vhr_out=vhtr, h_out=h_new, x_first_in=x_before_y) + x_before_y = .not. x_before_y do k=1,nz ; do j=jsd,jed ; do i=isd,ied - h_pre(i,j,k) = h_new(i,j,k)/G%areaT(i,j) !+ & -! max(0.0, 1.0e-13*h_vol(i,j,k) - G%areaT(i,j)*h_new(i,j,k))/G%areaT(i,j) + h_pre(i,j,k) = h_new(i,j,k)/G%areaT(i,j) enddo ; enddo ; enddo - + + + if(CS%debug) then + call hchksum(h_pre,"h_pre after advect_tracer",G%HI) + endif + call cpu_clock_begin(id_clock_ALE) - call ALE_main(G, GV, h_pre, CS%offline_CSp%u_preale, CS%offline_CSp%v_preale, CS%tv, & + call ALE_main_offline(G, GV, h_pre, CS%tv, & CS%tracer_Reg, CS%ALE_CSp, CS%offline_CSp%dt_offline) call cpu_clock_end(id_clock_ALE) + if(CS%debug) then + call hchksum(h_pre,"h_pre after ALE",G%HI) + endif + do k=1,nz ; do j=jsd,jed ; do i=isdB,iedB uhtr_sub(i,j,k) = uhtr(i,j,k) enddo ; enddo ; enddo @@ -1708,9 +1743,21 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) call pass_vector(uhtr_sub,vhtr_sub,G%Domain) call pass_var(h_pre, G%Domain) + + if(CS%debug) then + call uchksum(uhtr_sub,"uhtr_sub after adv iteration",G%HI) + call vchksum(vhtr_sub,"vhtr_sub after adv iteration",G%HI) + call hchksum(h_pre,"h_pre after adv iteration",G%HI) + endif - sum_u=sum(abs(uhtr)) - sum_v=sum(abs(vhtr)) + sum_u = 0.0 + do k=1,nz; do j=js,je ; do i=is-1,ie + sum_u = sum_u + abs(uhtr_sub(i,j,k)) + enddo; enddo; enddo + sum_v = 0.0 + do k=1,nz; do j=js-1,je; do i=is,ie + sum_v = sum_v + abs(vhtr_sub(i,j,k)) + enddo; enddo ; enddo call sum_across_PEs(sum_u) call sum_across_PEs(sum_v) @@ -1724,12 +1771,7 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) enddo ! Note here T/S are reset to the stored snap shot to ensure that the offline model ! densities, used in the neutral diffusion code don't drift too far from the online - ! model - CS%T = temp_old - CS%S = salt_old - call pass_var(CS%T,G%Domain) - call pass_var(CS%S,G%Domain) - + ! model call call_tracer_column_fns(h_pre, h_new, eatr*0.5, ebtr*0.5, & fluxes, CS%offline_CSp%dt_offline*0.5, G, GV, CS%tv, & CS%diabatic_CSp%optics, CS%tracer_flow_CSp, CS%debug, & @@ -1737,11 +1779,15 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) minimum_forcing_depth=0.001) call applyTracerBoundaryFluxesInOut(G, GV, zero_3dh, 0.5*dt_offline, fluxes, h_pre, & 0.8, 0.001) - + + if(CS%debug) then + call hchksum(h_pre,"h_pre after 2nd diabatic",G%HI) + endif + call cpu_clock_begin(id_clock_ALE) call ALE_offline_tracer_final( G, GV, h_pre, h_end, CS%tracer_Reg, CS%ALE_CSp) - + call cpu_clock_end(id_clock_ALE) endif - call pass_var(h_end, G%Domain) + ! Tracer diffusion Strang split between advection and diffusion call tracer_hordiff(h_end, CS%offline_CSp%dt_offline*0.5, CS%MEKE, CS%VarMix, G, GV, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv, CS%do_online, khdt_x*0.5, khdt_y*0.5) @@ -1759,10 +1805,14 @@ subroutine step_tracers(fluxes, state, Time_start, time_interval, CS) call disable_averaging(CS%diag) do i = is, ie ; do j = js, je ; do k=1,nz + CS%T(i,j,k) = temp_old(i,j,k) + CS%S(i,j,k) = salt_old(i,j,k) CS%h(i,j,k) = h_end(i,j,k) enddo ; enddo; enddo call pass_var(CS%h,G%Domain) + call pass_var(CS%T,G%Domain) + call pass_var(CS%S,G%Domain) diff --git a/src/tracer/MOM_offline_control.F90 b/src/tracer/MOM_offline_control.F90 index 4398b2de20..8e2f3b147f 100644 --- a/src/tracer/MOM_offline_control.F90 +++ b/src/tracer/MOM_offline_control.F90 @@ -209,9 +209,9 @@ subroutine transport_by_files(G, GV, CS, h_end, eatr, ebtr, uhtr, vhtr, khdt_x, timelevel=CS%ridx_mean,position=CENTER) !! Time-averaged fields - call read_data(CS%snap_file, 'temp_preadv', temp, domain=G%Domain%mpp_domain, & + call read_data(CS%snap_file, 'temp', temp, domain=G%Domain%mpp_domain, & timelevel=CS%ridx_mean,position=CENTER) - call read_data(CS%snap_file, 'salt_preadv', salt, domain=G%Domain%mpp_domain, & + call read_data(CS%snap_file, 'salt', salt, domain=G%Domain%mpp_domain, & timelevel=CS%ridx_mean,position=CENTER) !! Read snapshot fields (end of time interval timestamp) @@ -257,25 +257,6 @@ subroutine transport_by_files(G, GV, CS, h_end, eatr, ebtr, uhtr, vhtr, khdt_x, timelevel=CS%ridx_snap,position=center) endif -! if (do_ale) then -! CS%h_preale = GV%Angstrom -! CS%T_preale = 0.0 -! CS%S_preale = 0.0 -! CS%u_preale = 0.0 -! CS%v_preale = 0.0 -! call read_data(CS%preale_file, 'h_preale', CS%h_preale, domain=G%Domain%mpp_domain, & -! timelevel=CS%ridx_snap,position=CENTER) -! call read_data(CS%preale_file, 'T_preale', CS%T_preale, domain=G%Domain%mpp_domain, & -! timelevel=CS%ridx_mean,position=CENTER) -! call read_data(CS%preale_file, 'S_preale', CS%S_preale, domain=G%Domain%mpp_domain, & -! timelevel=CS%ridx_mean,position=CENTER) -! call read_data(CS%preale_file, 'u_preale', CS%u_preale, domain=G%Domain%mpp_domain, & -! timelevel=CS%ridx_mean,position=EAST) -! call read_data(CS%preale_file, 'v_preale', CS%v_preale, domain=G%Domain%mpp_domain, & -! timelevel=CS%ridx_mean,position=NORTH) -! -! endif - !! Make sure all halos have been updated ! Vector fields call pass_vector(uhtr, vhtr, G%Domain) @@ -291,13 +272,6 @@ subroutine transport_by_files(G, GV, CS, h_end, eatr, ebtr, uhtr, vhtr, khdt_x, if (do_ale) then call pass_var(fluxes%netMassOut,G%Domain) call pass_var(fluxes%netMassIn,G%Domain) -! -! call pass_vector(CS%u_preale,CS%v_preale,G%Domain) -! call pass_var(CS%h_preale, G%Domain) -! call pass_var(CS%T_preale, G%Domain) -! call pass_var(CS%S_preale, G%Domain) -! -! endif ! Update the read indices @@ -329,10 +303,6 @@ subroutine register_diags_offline_transport(Time, diag, CS) 'Meridional thickness fluxes remaining at end of timestep', 'kg') ! T-cell fields - CS%id_temp_preadv = register_diag_field('ocean_model', 'temp_preadv', diag%axesTL, Time, & - 'Temperature prior to advection', 'C') - CS%id_salt_preadv = register_diag_field('ocean_model', 'salt_preadv', diag%axesTL, Time, & - 'Salinity prior to advection', 'S') CS%id_hr = register_diag_field('ocean_model', 'hdiff', diag%axesTL, Time, & 'Difference between the stored and calculated layer thickness', 'm') CS%id_ear = register_diag_field('ocean_model', 'ear', diag%axesTL, Time, & @@ -606,232 +576,4 @@ subroutine limit_mass_flux_3d(G, GV, uh, vh, ea, eb, h_pre, max_off_cfl) end subroutine limit_mass_flux_3d - subroutine limit_mass_flux_ordered_3d(G, GV, uh, vh, ea, eb, h_in, h_end, max_off_cfl, z_first, x_before_y) - type(ocean_grid_type), pointer :: G - type(verticalGrid_type), pointer :: GV - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uh - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vh - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: ea - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: eb - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(in) :: h_in - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h_end - real, intent(in) :: max_off_cfl - logical, intent(in) :: z_first, x_before_y - - ! Local variables - integer :: i, j, k, m, is, ie, js, je, nz - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_budget ! Tracks how much thickness - ! remains for other fluxes - integer :: flux_order = -1 - ! In this subroutine, fluxes out of the box are scaled down if they deplete - ! the layer. Here the positive direction is defined as flux out of the box as opposed to the - ! typical, strictly upwind convention. Hence, uh(I-1) is multipled by negative one, - ! but uh(I) is not. This routine differs from limit_mass_flux_3d because in this case, - ! the ordering of direction matters. While this is more aggressive than the other routine which, - ! Scales fluxes if they would deplete the layer (independent of any convergence within an - ! iteration), this routine should still maintain a CFL less than 1 - ! Because horizontal transport must always be together (i.e. cannot do x->z->y), - ! four cases are considered) - ! 1: z -> x -> y - ! 2: z -> y -> x - ! 3: x -> y -> z - ! 4: y -> x -> z - - ! Set index-related variables for fields on T-grid - is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed ; nz = GV%ke - ! Copy layer thicknesses into a working array for this subroutine - do k=1,nz ; do j=js,je ; do i=is,ie - h_budget(i,j,k) = h_in(i,j,k)*G%areaT(i,j) - enddo ; enddo ; enddo - - ! Set the flux order (corresponding to one of the four cases described previously) - if (z_first .and. x_before_y) flux_order = 1 - if (z_first .and. (.not. x_before_y)) flux_order = 2 - if ((.not. z_first) .and. x_before_y) flux_order = 3 - if ((.not. z_first) .and. (.not. x_before_y)) flux_order = 4 - - select case (flux_order) - case (1) ! z -> x -> y - ! Check first to see if either the top or bottom flux would deplete the layer - !call flux_limiter_vertical(G, GV, h_budget, ea, eb, max_off_cfl) - call flux_limiter_u(G, GV, h_budget, uh, max_off_cfl) - call flux_limiter_v(G, GV, h_budget, vh, max_off_cfl) - case (2) ! z -> y -> x - !call flux_limiter_vertical(G, GV, h_budget, ea, eb, max_off_cfl) - call flux_limiter_v(G, GV, h_budget, vh, max_off_cfl) - call flux_limiter_u(G, GV, h_budget, uh, max_off_cfl) - case (3) ! x -> y -> z - call flux_limiter_u(G, GV, h_budget, uh, max_off_cfl) - call flux_limiter_v(G, GV, h_budget, vh, max_off_cfl) -! call flux_limiter_vertical(G, GV, h_budget, ea, eb, max_off_cfl) - case (4) ! y -> x -> z - call flux_limiter_v(G, GV, h_budget, vh, max_off_cfl) - call flux_limiter_u(G, GV, h_budget, uh, max_off_cfl) -! call flux_limiter_vertical(G, GV, h_budget, ea, eb, max_off_cfl) - case default - call MOM_error(FATAL, "Invalid choice of flux_order") - end select - - do k=1,nz ; do j=js,je ; do i=is,ie - h_end(i,j,k) = h_budget(i,j,k)/G%areaT(i,j) - if (h_end(i,j,k)<0.0 ) then - print *, "i,j,k,h,",i,j,k,h_end(i,j,k) - endif - enddo ; enddo ; enddo - - end subroutine limit_mass_flux_ordered_3d - - subroutine flux_limiter_vertical(G, GV, h, ea, eb, max_off_cfl) - type(ocean_grid_type), pointer :: G - type(verticalGrid_type), pointer :: GV - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: ea - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: eb - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h - real :: max_off_cfl - ! Limits how much the a layer can be depleted in the vertical direction - real, dimension(SZI_(G),SZK_(G)) :: ea2d, eb2d - real, dimension(SZI_(G),SZK_(G)) :: h2d, scale - real, dimension(SZI_(G),SZK_(G)+1) :: flux_interface, top_flux, bottom_flux - real :: total_out_flux, h_budget - integer :: i, j, k, m, is, ie, js, je, nz - - ! Set index-related variables for fields on T-grid - is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed ; nz = GV%ke - - do j=js,je - do k=1,nz ; do i=is,ie - ea2d(i,k) = ea(i,j,k) - eb2d(i,k) = eb(i,j,k) - h2d(i,k) = h(i,j,k) - enddo ; enddo; - - ! Calculate the fluxes through top and bottom faces of the cell - k=1 ! Top layer - do i=is,ie - flux_interface(i,k) = 0.0 - enddo - - do k=2, nz-1 ; do i=is,ie - top_flux(i,k) = ea2d(i,k)-eb2d(i,k-1) - bottom_flux(i,k) = ea2d(i,k+1)-eb2d(i,k) - enddo ; enddo - k=nz ! Bottom layer - do i=is,ie - top_flux(i,k) = ea2d(i,k)-eb2d(i,k-1) - bottom_flux(i,k) = -eb2d(i,k) - enddo - - ! Convert fluxes which are in units of thickness to units of volume - do k=1,nz ; do i=is,ie - top_flux(i,k) = top_flux(i,k)*G%areaT(i,j) - bottom_flux(i,k) = bottom_flux(i,k)*G%areaT(i,j) - enddo; enddo - - do k=1,nz ; do i=is,ie - - enddo; enddo - - enddo - end subroutine flux_limiter_vertical - - subroutine flux_limiter_u(G, GV, h, uh, max_off_cfl) - type(ocean_grid_type), pointer :: G - type(verticalGrid_type), pointer :: GV - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uh - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h - real :: max_off_cfl - ! Limits how much the a layer can be depleted in the vertical direction - real, dimension(SZIB_(G),SZK_(G)) :: uh2d - real :: hup, hlos, min_h, h_remain - integer :: i, j, k, m, is, ie, js, je, nz - - min_h= 0.1*GV%Angstrom - ! Set index-related variables for fields on T-grid - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - - do j=js-1,je - do k=1,nz ; do i=is-2,ie - uh2d(I,k) = uh(I,j,k) - enddo ; enddo; - - do k=1,nz ; do i=is-1,ie - if(uh2d(I,k)<0.0) then - hup = h(i+1,j,k) - min_h*G%areaT(i+1,j) - hlos = max(0.0, uh2d(I+1,k)) - if (( ((hup-hlos)+uh2d(I,k))<0.0) .and. & - ((0.5*hup + uh2d(I,k))<0.0)) then - uh2d(I,k) = min(-0.5*hup,-hup+hlos,0.0) - endif - else - hup = h(i,j,k) - G%areaT(i,j)*min_h - hlos = max(0.0,-uh2d(I-1,k)) - if ((((hup-hlos)-uh2d(I,k))<0.0) .and. & - ((0.5*hup-uh2d(I,k))<0.0)) then - uh2d(I,k) = max(0.5*hup,hup-hlos,0.0) - endif - endif - enddo ; enddo - - do k=1,nz - do i=is-1,ie - uh(I,j,k) = uh2d(I,k) - enddo - do i=is,ie - h(i,j,k) = h(i,j,k) - (uh2d(I,k) - uh2d(I-1,k)) - enddo - enddo - - enddo - - end subroutine flux_limiter_u - - subroutine flux_limiter_v(G, GV, h, vh, max_off_cfl) - type(ocean_grid_type), pointer :: G - type(verticalGrid_type), pointer :: GV - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vh - real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h - real, intent(in) :: max_off_cfl - ! Limits how much the a layer can be depleted in the vertical direction - real, dimension(SZJB_(G),SZK_(G)) :: vh2d - real, dimension(SZJ_(G),SZK_(G)) :: h2d - real :: hup, hlos, min_h, h_remain - integer :: i, j, k, m, is, ie, js, je, nz - - ! Set index-related variables for fields on T-grid - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - min_h= 0.1*GV%Angstrom - do i=is-1,ie - do k=1,nz ; do j=js-2,je - vh2d(J,k) = vh(i,J,k) - enddo ; enddo; - do k=1,nz ; do j=js-1,je - if(vh2d(J,k)<0.0) then - hup = h(i,j+1,k)-G%areaT(i,j+1)*min_h - hlos = max(0.0,vh2d(J+1,k)) - if ((((hup-hlos)+vh2d(J,k))<0.0) .and. & - ((0.5*hup+vh2d(J,k))<0.0)) then - vh2d(J,k) = min(-0.5*hup,-hup+hlos,0.0) - endif - else - hup = h(i,j,k) - G%areaT(i,j)*min_h - hlos = max(0.0,-vh2d(J-1,k)) - if ((((hup-hlos)-vh2d(J,k))<0.0) .and. & - ((0.5*hup - vh2d(J,k))<0.0)) then - vh2d(J,k) = max(0.5*hup,hup-hlos,0.0) - endif - endif - enddo ; enddo - - do k=1,nz - do j=js-1,je - vh(i,J,k) = vh2d(J,k) - enddo - do j=js,je - h(i,j,k) = h(i,j,k) - (vh2d(J,k) - vh2d(J-1,k)) - enddo - enddo - enddo - end subroutine flux_limiter_v - - end module MOM_offline_transport