From 024157fadf9edb20a6cfe0312c2b0512b00e7829 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 27 Oct 2016 19:04:28 -0400 Subject: [PATCH 1/3] Eliminates the need to use ZSTAR_RIGID_SURFACE_THRESHOLD This commit eliminates the need to specify ZSTAR_RIGID_SURFACE_THRESHOLD when using ALE/ZSTAR with an ice shelf. --- src/ALE/MOM_ALE.F90 | 25 ++++++-- src/ALE/MOM_regridding.F90 | 29 ++++++--- src/core/MOM.F90 | 63 +++++++++++++++++-- .../MOM_state_initialization.F90 | 44 ++++++++++++- 4 files changed, 140 insertions(+), 21 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index fa350742b6..88fccf1857 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -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) @@ -381,13 +381,17 @@ 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 nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec + ice_shelf = associated(frac_shelf_h) + if (CS%show_call_tree) call callTree_enter("ALE_main(), MOM_ALE.F90") if (present(dt)) then @@ -397,7 +401,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. ) @@ -554,7 +562,7 @@ 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 @@ -562,20 +570,25 @@ subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug ) 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. diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 3b39e753a6..9daca62674 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -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 @@ -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 ) @@ -571,7 +579,7 @@ 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 @@ -579,16 +587,18 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface ) 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 nz = GV%ke minThickness = CS%min_thickness + ice_shelf = associated(frac_shelf_h) !$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h) & !$OMP private(nominalDepth,totalThickness,minThickness, & @@ -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 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 @@ -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 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7de2ce4c1b..e2bd1ff5a7 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -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 @@ -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 @@ -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 @@ -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 )" @@ -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 @@ -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 @@ -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 @@ -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 "// & @@ -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. @@ -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) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 8e4f014830..503545cd9d 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1608,9 +1608,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) type(directories), intent(in) :: dirs character(len=200) :: filename ! The name of an input file containing temperature - ! and salinity in z-space. + ! and salinity in z-space; also used for ice shelf area. character(len=200) :: inputdir ! The directory where NetCDF input files are. - character(len=200) :: mesg + character(len=200) :: mesg, area_varname, ice_shelf_file type(EOS_type), pointer :: eos => NULL() @@ -1630,6 +1630,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) integer :: kd, inconsistent real :: PI_180 ! for conversion from degrees to radians + real, dimension(:,:), pointer :: shelf_area real :: min_depth real :: dilate real :: missing_value_temp, missing_value_salt @@ -1657,12 +1658,15 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) ! Local variables for ALE remapping real, dimension(:), allocatable :: hTarget + real, dimension(:,:), allocatable :: area_shelf_h + real, dimension(:,:), allocatable, target :: frac_shelf_h real, dimension(:,:,:), allocatable :: tmpT1dIn, tmpS1dIn, h1, tmp_mask_in real :: zTopOfCell, zBottomOfCell type(regridding_CS) :: regridCS ! Regridding parameters and work arrays type(remapping_CS) :: remapCS ! Remapping parameters and work arrays logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg + logical :: use_ice_shelf character(len=10) :: remappingScheme real :: tempAvg, saltAvg integer :: nPoints, ans @@ -1758,6 +1762,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) kd = size(z_in,1) allocate(rho_z(isd:ied,jsd:jed,kd)) + allocate(area_shelf_h(isd:ied,jsd:jed)) + allocate(frac_shelf_h(isd:ied,jsd:jed)) press(:)=tv%p_ref @@ -1772,6 +1778,33 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) call pass_var(mask_z,G%Domain) call pass_var(rho_z,G%Domain) + ! This is needed for building an ALE grid under ice shelves + call get_param(PF, mod, "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.) + if (use_ice_shelf) then + call get_param(PF, mod, "ICE_THICKNESS_FILE", ice_shelf_file, & + "The file from which the ice bathymetry and area are read.", & + fail_if_missing=.true.) + call log_param(PF, mod, "INPUTDIR/THICKNESS_FILE", filename) + call get_param(PF, mod, "ICE_AREA_VARNAME", area_varname, & + "The name of the area variable in ICE_THICKNESS_FILE.", & + fail_if_missing=.true.) + if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & + "MOM_temp_salt_initialize_from_Z: Unable to open "//trim(filename)) + + 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 + + endif + ! Done with horizontal interpolation. ! Now remap to model coordinates if (useALEremapping) then @@ -1853,7 +1886,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) call pass_var(h, G%Domain) ! Regridding might eventually use spatial information and call pass_var(tv%T, G%Domain) ! thus needs to be up to date in the halo regions even though call pass_var(tv%S, G%Domain) ! ALE_build_grid() only updates h on the computational domain. - call ALE_build_grid( G, GV, regridCS, remapCS, h, tv, .true. ) + + if (use_ice_shelf) then + call ALE_build_grid( G, GV, regridCS, remapCS, h, tv, .true., shelf_area) + else + call ALE_build_grid( G, GV, regridCS, remapCS, h, tv, .true. ) + endif endif call ALE_remap_scalar( remapCS, G, GV, nz, h1, tmpT1dIn, h, tv%T, all_cells=remap_full_column, old_remap=remap_old_alg ) call ALE_remap_scalar( remapCS, G, GV, nz, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, old_remap=remap_old_alg ) From 136faceb695d06279af8efb82a290679d436e932 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 28 Oct 2016 08:57:39 -0400 Subject: [PATCH 2/3] Deleted the get_param call for ZSTAR_RIGID_SURFACE_THRESHOLD. --- src/ALE/MOM_ALE.F90 | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 88fccf1857..b31c04d212 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -1042,7 +1042,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. @@ -1077,15 +1076,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"//& From 52efd2cb774bb2e53a12f825b17ad3678e77c97b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 28 Oct 2016 11:11:06 -0400 Subject: [PATCH 3/3] Adds ice_shelf = false as the default value --- src/ALE/MOM_ALE.F90 | 5 +++-- src/ALE/MOM_regridding.F90 | 8 ++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index b31c04d212..486061d39b 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -386,11 +386,12 @@ subroutine ALE_main( G, GV, h, u, v, tv, Reg, CS, dt, frac_shelf_h) 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 + logical :: ice_shelf = .false. nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec - ice_shelf = associated(frac_shelf_h) + 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") diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 9daca62674..1887beae35 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -594,15 +594,15 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) real :: nominalDepth, totalThickness, dh real, dimension(SZK_(GV)+1) :: zOld, zNew real :: minThickness - logical :: ice_shelf + logical :: ice_shelf = .false. nz = GV%ke minThickness = CS%min_thickness - ice_shelf = associated(frac_shelf_h) + 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