Skip to content

Commit

Permalink
Answers do not change across PE counts
Browse files Browse the repository at this point in the history
  • Loading branch information
Andrew Shao committed Oct 6, 2016
1 parent b21e331 commit 724dc3e
Show file tree
Hide file tree
Showing 3 changed files with 132 additions and 284 deletions.
56 changes: 56 additions & 0 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
98 changes: 74 additions & 24 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)


Expand Down Expand Up @@ -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
Expand All @@ -1662,22 +1676,34 @@ 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
do k=1,nz ; do j=jsdB,jedB ; do i=isd,ied
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, &
evap_CFL_limit=0.8, &
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

Expand All @@ -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
Expand All @@ -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)
Expand All @@ -1724,24 +1771,23 @@ 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, &
evap_CFL_limit=0.8, &
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)
Expand All @@ -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)



Expand Down
Loading

0 comments on commit 724dc3e

Please sign in to comment.