Skip to content

Eliminates the need to use ZSTAR_RIGID_SURFACE_THRESHOLD #378

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
adcroft merged 3 commits into from
Nov 1, 2016
Merged

Eliminates the need to use ZSTAR_RIGID_SURFACE_THRESHOLD #378

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
024157f
Eliminates the need to use ZSTAR_RIGID_SURFACE_THRESHOLD
gustavo-marques Oct 27, 2016
136face
Deleted the get_param call for ZSTAR_RIGID_SURFACE_THRESHOLD.
gustavo-marques Oct 28, 2016
52efd2c
Adds ice_shelf = false as the default value
gustavo-marques Oct 28, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 20 additions & 16 deletions src/ALE/MOM_ALE.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -371,7 +371,7 @@ end subroutine ALE_end
!! 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( G, GV, h, u, v, tv, Reg, CS, dt)
subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt, frac_shelf_h)
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)
Expand All @@ -381,13 +381,18 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt)
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()
real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
! 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
logical :: ice_shelf = .false.

nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec

if (present(frac_shelf_h)) ice_shelf = .true.
!ice_shelf = associated(frac_shelf_h)

if (CS%show_call_tree) call callTree_enter("ALE_main(), MOM_ALE.F90")

if (present(dt)) then
Expand All @@ -397,7 +402,11 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt)

! 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 )
if (ice_shelf) then
call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, frac_shelf_h)
else
call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid)
endif

call check_grid( G, GV, h, 0. )

Expand Down Expand Up @@ -554,28 +563,33 @@ subroutine check_grid( G, GV, h, threshold )
end subroutine check_grid

!> Generates new grid
subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug )
subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h )
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(regridding_CS), intent(in) :: regridCS !< Regridding parameters and options
type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options
type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variable structure
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the last time step (m or Pa)
logical, optional, intent(in) :: debug !< If true, show the call tree

real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
! Local variables
integer :: nk, i, j, k
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 ! The new grid thicknesses
logical :: show_call_tree
logical :: show_call_tree, use_ice_shelf

show_call_tree = .false.
if (present(frac_shelf_h)) use_ice_shelf = .true.
if (present(debug)) show_call_tree = debug
if (show_call_tree) call callTree_enter("ALE_build_grid(), MOM_ALE.F90")

! 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( remapCS, regridCS, G, GV, h, tv, h_new, dzRegrid )
if (use_ice_shelf) then
call regridding_main( remapCS, regridCS, G, GV, h, tv, h_new, dzRegrid, frac_shelf_h )
else
call regridding_main( remapCS, regridCS, G, GV, h, tv, h_new, dzRegrid )
endif

! Override old grid with new one. The new grid 'h_new' is built in
! one of the 'build_...' routines above.
Expand Down Expand Up @@ -1029,7 +1043,6 @@ subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz )
real :: filt_len, strat_tol, index_scale
real :: tmpReal, compress_fraction
real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
real :: height_of_rigid_surface
integer :: nz_fixed_sfc
real :: rho_target(GV%ke+1) ! Target density used in HYBRID mode
real, dimension(size(dz)) :: h_max ! Maximum layer thicknesses, in m.
Expand Down Expand Up @@ -1064,15 +1077,6 @@ subroutine ALE_initRegridding(GV, max_depth, param_file, mod, regridCS, dz )
call initialize_regridding( GV%ke, coordMode, interpScheme, regridCS, &
compressibility_fraction=compress_fraction )

if (coordMode(1:2) == 'Z*') then
call get_param(param_file, mod, "ZSTAR_RIGID_SURFACE_THRESHOLD", height_of_rigid_surface, &
"A threshold height used to detect the presence of a rigid-surface\n"//&
"depressing the upper-surface of the model, such as an ice-shelf.\n"//&
"This is a temporary work around for initialization under an ice-shelf.", &
units='m', default=-1.E30)
call set_regrid_params( regridCS, height_of_rigid_surface=height_of_rigid_surface*GV%m_to_H)
endif

call get_param(param_file, mod, "ALE_COORDINATE_CONFIG", string, &
"Determines how to specify the coordinate\n"//&
"resolution. Valid options are:\n"//&
Expand Down
33 changes: 23 additions & 10 deletions src/ALE/MOM_regridding.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -266,7 +266,7 @@ end subroutine end_regridding
!------------------------------------------------------------------------------
! Dispatching regridding routine: regridding & remapping
!------------------------------------------------------------------------------
subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface )
subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h)
!------------------------------------------------------------------------------
! This routine takes care of (1) building a new grid and (2) remapping between
! the old grid and the new grid. The creation of the new grid can be based
Expand All @@ -293,14 +293,22 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface )
type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...)
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h_new !< New 3D grid consistent with target coordinate
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in position of each interface
real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage
! Local variables
real :: trickGnuCompiler
logical :: use_ice_shelf

if (present(frac_shelf_h)) use_ice_shelf = .true.

select case ( CS%regridding_scheme )

case ( REGRIDDING_ZSTAR )
call build_zstar_grid( CS, G, GV, h, dzInterface )
call calc_h_new_by_dz(G, GV, h, dzInterface, h_new)
if (use_ice_shelf) then
call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h )
else
call build_zstar_grid( CS, G, GV, h, dzInterface )
endif
call calc_h_new_by_dz(G, GV, h, dzInterface, h_new)

case ( REGRIDDING_SIGMA )
call build_sigma_grid( CS, G, GV, h, dzInterface )
Expand Down Expand Up @@ -571,28 +579,30 @@ end subroutine filtered_grid_motion
!> Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004).
!! z* is defined as
!! z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H .
subroutine build_zstar_grid( CS, G, GV, h, dzInterface )
subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)

! Arguments
type(regridding_CS), intent(in) :: CS !< Regridding control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses, in H
real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth in H.

real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage.
! Local variables
integer :: i, j, k
integer :: nz
real :: nominalDepth, totalThickness, dh
real, dimension(SZK_(GV)+1) :: zOld, zNew
real :: minThickness
logical :: ice_shelf = .false.

nz = GV%ke
minThickness = CS%min_thickness
if (present(frac_shelf_h)) ice_shelf = .true.

!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h) &
!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h) &
!$OMP private(nominalDepth,totalThickness,minThickness, &
!$OMP zNew,dh,zOld)
!$OMP zNew,dh,zOld,ice_shelf)
do j = G%jsc-1,G%jec+1
do i = G%isc-1,G%iec+1

Expand All @@ -615,10 +625,14 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface )
zOld(k) = zOld(k+1) + h(i,j,k)
enddo

if (totalThickness-nominalDepth<CS%height_of_rigid_surface) then
call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew, &
if (ice_shelf) then
if (frac_shelf_h(i,j) > 0.) then ! under ice shelf
call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew, &
z_rigid_top = totalThickness-nominalDepth, &
eta_orig = zOld(1))
else
call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew)
endif
else
call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew)
endif
Expand Down Expand Up @@ -2782,7 +2796,6 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
CS%halocline_strat_tol = halocline_strat_tol
endif
if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e
if (present(height_of_rigid_surface)) CS%height_of_rigid_surface = height_of_rigid_surface

end subroutine set_regrid_params

Expand Down
63 changes: 59 additions & 4 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ module MOM
use MOM_forcing_type, only : MOM_forcing_chksum
use MOM_get_input, only : Get_MOM_Input, directories
use MOM_io, only : MOM_io_init, vardesc, var_desc
use MOM_io, only : slasher, file_exists, read_data
use MOM_obsolete_params, only : find_obsolete_params
use MOM_restart, only : register_restart_field, query_initialized, save_restart
use MOM_restart, only : restart_init, MOM_restart_CS
Expand Down Expand Up @@ -494,7 +495,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
type(time_type) :: Time_local
logical :: showCallTree
logical :: do_pass_kd_kv_turb ! This is used for a group halo pass.

logical :: use_ice_shelf ! Needed for ALE
G => CS%G ; GV => CS%GV
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
Expand All @@ -508,6 +509,8 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
showCallTree = callTree_showQuery()
if (showCallTree) call callTree_enter("step_MOM(), MOM.F90")

if (associated(fluxes%frac_shelf_h)) use_ice_shelf = .true.

! First determine the time step that is consistent with this call.
! It is anticipated that the time step will almost always coincide
! with dt. In addition, ntstep is determined, subject to the constraint
Expand Down Expand Up @@ -741,7 +744,14 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
call check_redundant("Pre-ALE 1 ", u, v, G)
endif
call cpu_clock_begin(id_clock_ALE)
call ALE_main(G, GV, h, u, v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, dtdia)
if (use_ice_shelf) then

call ALE_main(G, GV, h, u, v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, dtdia, &
fluxes%frac_shelf_h)
else
call ALE_main(G, GV, h, u, v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, dtdia)
endif

call cpu_clock_end(id_clock_ALE)
endif ! endif for the block "if ( CS%use_ALE_algorithm )"

Expand Down Expand Up @@ -1084,7 +1094,13 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
call check_redundant("Pre-ALE ", u, v, G)
endif
call cpu_clock_begin(id_clock_ALE)
call ALE_main(G, GV, h, u, v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%dt_trans)
if (use_ice_shelf) then

call ALE_main(G, GV, h, u, v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%dt_trans, &
fluxes%frac_shelf_h)
else
call ALE_main(G, GV, h, u, v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%dt_trans)
endif
call cpu_clock_end(id_clock_ALE)
endif

Expand Down Expand Up @@ -1930,6 +1946,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo

real, allocatable, dimension(:,:,:) :: e ! interface heights (meter)
real, allocatable, dimension(:,:) :: eta ! free surface height (m) or bottom press (Pa)
real, allocatable, dimension(:,:) :: area_shelf_h ! area occupied by ice shelf
real, dimension(:,:), allocatable, target :: frac_shelf_h ! fraction of total area occupied by ice shelf
real, dimension(:,:), pointer :: shelf_area
type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL()

real :: default_val ! default value for a parameter
Expand All @@ -1941,6 +1960,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo
logical :: save_IC ! If true, save the initial conditions.
logical :: do_unit_tests ! If true, call unit tests.
logical :: test_grid_copy = .false.
logical :: use_ice_shelf ! Needed for ALE
logical :: global_indexing ! If true use global horizontal index values instead
! of having the data domain on each processor start at 1.
logical :: bathy_at_vel ! If true, also define bathymetric fields at the
Expand All @@ -1953,6 +1973,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo

type(time_type) :: Start_time
type(ocean_internal_state) :: MOM_internal_state
character(len=200) :: area_varname, ice_shelf_file, inputdir, filename

if (associated(CS)) then
call MOM_error(WARNING, "initialize_MOM called with an associated "// &
Expand Down Expand Up @@ -2223,6 +2244,19 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo
"initialize_MOM: A bulk mixed layer can only be used with T & S as "//&
"state variables. Add USE_EOS = True to MOM_input.")

call get_param(param_file, 'MOM', "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)
if (use_ice_shelf) then
inputdir = "." ; call get_param(param_file, 'MOM', "INPUTDIR", inputdir)
inputdir = slasher(inputdir)
call get_param(param_file, 'MOM', "ICE_THICKNESS_FILE", ice_shelf_file, &
"The file from which the ice bathymetry and area are read.", &
fail_if_missing=.true.)
call get_param(param_file, 'MOM', "ICE_AREA_VARNAME", area_varname, &
"The name of the area variable in ICE_THICKNESS_FILE.", &
fail_if_missing=.true.)
endif


call callTree_waypoint("MOM parameters read (initialize_MOM)")

! Set up the model domain and grids.
Expand Down Expand Up @@ -2464,7 +2498,28 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo
call callTree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
call adjustGridForIntegrity(CS%ALE_CSp, G, GV, CS%h )
call callTree_waypoint("Calling ALE_main() to remap initial conditions (initialize_MOM)")
call ALE_main( G, GV, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp )
if (use_ice_shelf) then
filename = trim(inputdir)//trim(ice_shelf_file)
if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
"MOM: Unable to open "//trim(filename))

allocate(area_shelf_h(isd:ied,jsd:jed))
allocate(frac_shelf_h(isd:ied,jsd:jed))
call read_data(filename,trim(area_varname),area_shelf_h,domain=G%Domain%mpp_domain)
! initialize frac_shelf_h with zeros (open water everywhere)
frac_shelf_h(:,:) = 0.0
! compute fractional ice shelf coverage of h
do j=jsd,jed ; do i=isd,ied
if (G%areaT(i,j) > 0.0) &
frac_shelf_h(i,j) = area_shelf_h(i,j) / G%areaT(i,j)
enddo ; enddo
! pass to the pointer
shelf_area => frac_shelf_h
call ALE_main(G, GV, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, &
frac_shelf_h = shelf_area)
else
call ALE_main( G, GV, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp )
endif
call cpu_clock_begin(id_clock_pass_init)
call do_group_pass(CS%pass_uv_T_S_h, G%Domain)
call cpu_clock_end(id_clock_pass_init)
Expand Down
Loading