diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 0732e1b0b8..a97fdd7dce 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -705,7 +705,7 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c if (present(old_remap)) use_remapping_core_w = old_remap n_points = nk_src -!$OMP parallel default(none) shared(CS,G,h_src,s_src,h_dst,s_dst & +!$OMP parallel default(none) shared(CS,G,GV,h_src,s_src,h_dst,s_dst & !$OMP ,ignore_vanished_layers, use_remapping_core_w, nk_src,dx ) & !$OMP firstprivate(n_points) !$OMP do @@ -917,6 +917,7 @@ 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. @@ -951,6 +952,15 @@ 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"//& diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 0c8f98eb4e..b59cebd277 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -133,6 +133,13 @@ module MOM_regridding !! If false, integrate from the bottom upward, as does the rest of the model. logical :: integrate_downward_for_e = .true. + + !> The height to use as a threshold for detection of a "rigid lid" such as an ice shelf + !! base. If the free-surface is depressed below this height the z* coordinate generator + !! assumes the depression is due to a rigid surface. This is a kludge until we provide + !! the position of such a surface to the ALE code. + real :: height_of_rigid_surface = -1.E30 + end type ! The following routines are visible to the outside world @@ -337,7 +344,7 @@ subroutine calc_h_new_by_dz(G, GV, h, dzInterface, h_new) ! Local variables integer :: i, j, k -!$OMP parallel do default(none) shared(G,h,dzInterface,h_new) +!$OMP parallel do default(none) shared(G,GV,h,dzInterface,h_new) do j = G%jsc-1,G%jec+1 do i = G%isc-1,G%iec+1 if (G%mask2dT(i,j)>0.) then @@ -362,7 +369,7 @@ subroutine check_remapping_grid( G, GV, h, dzInterface, msg ) ! Local variables integer :: i, j -!$OMP parallel do default(none) shared(G,h,dzInterface,msg) +!$OMP parallel do default(none) shared(G,GV,h,dzInterface,msg) do j = G%jsc-1,G%jec+1 do i = G%isc-1,G%iec+1 if (G%mask2dT(i,j)>0.) call check_grid_column( GV%ke, G%bathyT(i,j), h(i,j,:), dzInterface(i,j,:), msg ) @@ -602,13 +609,19 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface ) totalThickness = totalThickness + h(i,j,k) end do - call build_zstar_column(CS, nz, nominalDepth, totalThickness, zNew) - zOld(nz+1) = - nominalDepth do k = nz,1,-1 zOld(k) = zOld(k+1) + h(i,j,k) enddo + if (totalThickness-nominalDepth Builds a z* coordinate with a minimum thickness -subroutine build_zstar_column( CS, nz, depth, total_thickness, zInterface) +subroutine build_zstar_column(CS, nz, depth, total_thickness, zInterface, z_rigid_top, eta_orig) type(regridding_CS), intent(in) :: CS !< Regridding control structure integer, intent(in) :: nz !< Number of levels real, intent(in) :: depth !< Depth of ocean bottom (positive in m) real, intent(in) :: total_thickness !< Column thickness (positive in m) real, dimension(nz+1), intent(inout) :: zInterface !< Absolute positions of interfaces + real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (negative in m) + real, optional, intent(in) :: eta_orig !< The actual original height of the top (m) ! Local variables - real :: eta, stretching, dh, min_thickness + real :: eta, stretching, dh, min_thickness, z0_top, z_star integer :: k + logical :: new_zstar_def + new_zstar_def = .false. min_thickness = min( CS%min_thickness, total_thickness/real(nz) ) + z0_top = 0. + if (present(z_rigid_top)) then + z0_top = z_rigid_top + new_zstar_def = .true. + endif - ! Position of free-surface + ! Position of free-surface (or the rigid top, for which eta ~ z0_top) eta = total_thickness - depth + if (present(eta_orig)) eta = eta_orig + + ! Conventional z* coordinate: + ! z* = (z-eta) / stretching where stretching = (H+eta)/H + ! z = eta + stretching * z* + ! The above gives z*(z=eta) = 0, z*(z=-H) = -H. + ! With a rigid top boundary at eta = z0_top then + ! z* = z0 + (z-eta) / stretching where stretching = (H+eta)/(H+z0) + ! z = eta + stretching * (z*-z0) * stretching + stretching = total_thickness / ( depth + z0_top ) + + if (new_zstar_def) then + ! z_star is the notional z* coordinate in absence of upper/lower topography + z_star = 0. ! z*=0 at the free-surface + zInterface(1) = eta ! The actual position of the top of the column + do k = 2,nz + z_star = z_star - CS%coordinateResolution(k-1) + ! This ensures that z is below a rigid upper surface (ice shelf bottom) + zInterface(k) = min( eta + stretching * ( z_star - z0_top ), z0_top ) + ! This ensures that the layer in inflated + zInterface(k) = min( zInterface(k), zInterface(k-1) - min_thickness ) + ! This ensures that z is above or at the topography + zInterface(k) = max( zInterface(k), -depth + real(nz+1-k) * min_thickness ) + enddo + zInterface(nz+1) = -depth - ! z* = (z-eta) / stretching where stretching = (H+eta)/H - ! z = eta + stretching * z* - stretching = total_thickness / depth - - ! Integrate down from the top for a notional new grid, ignoring topography - zInterface(1) = eta - do k = 1,nz - dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing - zInterface(k+1) = zInterface(k) - dh - enddo + else + ! Integrate down from the top for a notional new grid, ignoring topography + ! The starting position is offset by z0_top which, if z0_top<0, will place + ! interfaces above the rigid boundary. + zInterface(1) = eta + do k = 1,nz + dh = stretching * CS%coordinateResolution(k) ! Notional grid spacing + zInterface(k+1) = zInterface(k) - dh + enddo - ! Integrating up from the bottom adjusting interface position to accommodate - ! inflating layers without disturbing the interface above - zInterface(nz+1) = -depth - do k = nz,1,-1 - if ( zInterface(k) < (zInterface(k+1) + min_thickness) ) then - zInterface(k) = zInterface(k+1) + min_thickness - endif - enddo + ! Integrating up from the bottom adjusting interface position to accommodate + ! inflating layers without disturbing the interface above + zInterface(nz+1) = -depth + do k = nz,1,-1 + if ( zInterface(k) < (zInterface(k+1) + min_thickness) ) then + zInterface(k) = zInterface(k+1) + min_thickness + endif + enddo + endif end subroutine build_zstar_column @@ -1699,7 +1746,7 @@ subroutine adjust_interface_motion( nk, min_thickness, h_old, dz_int ) 'implied h<0 is larger than roundoff!') endif enddo - do k = nk,1,-1 + do k = nk,2,-1 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) ) if (h_new0.) &" is G was valid on BT domain if (CS%bathyT(i+1,j)+CS%bathyT(i,j)>0.) & Datu(I,j) = 2.0*CS%dy_Cu(I,j) * GV%m_to_H * & @@ -3549,6 +3550,7 @@ subroutine find_face_areas(Datu, Datv, G, GV, CS, MS, eta, halo, add_max) enddo ; enddo !$OMP do do J=js-1-hs,je+hs ; do i=is-hs,ie+hs + Datv(i, J) = 0.0 !Would be "if (G%mask2dCv(i,J)>0.) &" is G was valid on BT domain if (CS%bathyT(i,j+1)+CS%bathyT(i,j)>0.) & Datv(i,J) = 2.0*CS%dx_Cv(i,J) * GV%m_to_H * & diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0f9568cd2c..8121d5af6a 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -77,6 +77,7 @@ module MOM_state_initialization use MOM_ALE, only : ALE_remap_scalar, ALE_build_grid use MOM_regridding, only : regridding_CS, set_regrid_params use MOM_remapping, only : remapping_CS, initialize_remapping +use MOM_remapping, only : remapping_core_h use MOM_tracer_initialization_from_Z, only : horiz_interp_and_extrap_tracer implicit none ; private @@ -846,7 +847,7 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(ALE_CS), pointer :: ALE_CSp !< ALE control structure - type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness (H units, m or Pa) ! Local variables character(len=200) :: mod = "trim_for_ice" @@ -856,6 +857,8 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path real :: scale_factor, min_thickness integer :: i, j, k + logical :: use_remapping + type(remapping_CS), pointer :: remap_CS => NULL() call get_param(PF, mod, "SURFACE_PRESSURE_FILE", p_surf_file, & "The initial condition file for the surface height.", & @@ -875,6 +878,13 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) if (scale_factor /= 1.) p_surf(:,:) = scale_factor * p_surf(:,:) call get_param(PF, mod, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', & units='m', default=1.e-3) + call get_param(PF, mod, "TRIMMING_USES_REMAPPING", use_remapping, & + 'When trimming the column, also remap T and S.', & + default=.false.) + if (use_remapping) then + allocate(remap_CS) + call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true.) + endif ! Find edge values of T and S used in reconstructions if ( associated(ALE_CSp) ) then ! This should only be associated if we are in ALE mode @@ -893,30 +903,34 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) do j=G%jsc,G%jec ; do i=G%isc,G%iec call cut_off_column_top(GV%ke, tv, GV%Rho0, G%G_earth, G%bathyT(i,j), min_thickness, & - tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), p_surf(i,j), h(i,j,:)) + tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), & + p_surf(i,j), h(i,j,:), remap_CS) enddo ; enddo end subroutine trim_for_ice !> Adjust the layer thicknesses by cutting away the top at the depth where the hydrostatic !! pressure matches p_surf -subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, T, T_t, T_b, S, S_t, S_b, p_surf, h) +subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, & + T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS) integer, intent(in) :: nk !< Number of layers type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure real, intent(in) :: Rho0 !< Reference density (kg/m3) real, intent(in) :: G_earth !< Gravitational acceleration (m/s2) real, intent(in) :: depth !< Depth of ocean column (m) real, intent(in) :: min_thickness !< Smallest thickness allowed (m) - real, dimension(nk), intent(in) :: T !< - real, dimension(nk), intent(in) :: T_t !< - real, dimension(nk), intent(in) :: T_b !< - real, dimension(nk), intent(in) :: S !< - real, dimension(nk), intent(in) :: S_t !< - real, dimension(nk), intent(in) :: S_b !< + real, dimension(nk), intent(inout) :: T !< Layer mean temperature + real, dimension(nk), intent(in) :: T_t !< Temperature at top of layer + real, dimension(nk), intent(in) :: T_b !< Temperature at bottom of layer + real, dimension(nk), intent(inout) :: S !< Layer mean salinity + real, dimension(nk), intent(in) :: S_t !< Salinity at top of layer + real, dimension(nk), intent(in) :: S_b !< Salinity at bottom of layer real, intent(in) :: p_surf !< Imposed pressure on ocean at surface (Pa) real, dimension(nk), intent(inout) :: h !< Layer thickness (H units, m or Pa) + type(remapping_CS), pointer :: remap_CS ! Remapping structure for remapping T and S, if associated ! Local variables real, dimension(nk+1) :: e ! Top and bottom edge values for reconstructions + real, dimension(nk) :: h0, S0, T0, h1, S1, T1 real :: P_t, P_b, z_out, e_top integer :: k @@ -924,6 +938,7 @@ subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, T, T_ e(nk+1) = -depth do k=nk,1,-1 e(K) = e(K+1) + h(k) + h0(k) = h(nk+1-k) ! Keep a copy to use in remapping enddo P_t = 0. @@ -958,6 +973,22 @@ subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, T, T_ enddo endif + ! Now we need to remap but remapping assumes the surface is at the + ! same place in the two columns so we turn the column upside down. + if (associated(remap_CS)) then + do k=1,nk + S0(k) = S(nk+1-k) + T0(k) = T(nk+1-k) + h1(k) = h(nk+1-k) + enddo + call remapping_core_h(nk, h0, T0, nk, h1, T1, remap_CS ) + call remapping_core_h(nk, h0, S0, nk, h1, S1, remap_CS ) + do k=1,nk + S(k) = S1(nk+1-k) + T(k) = T1(nk+1-k) + enddo + endif + end subroutine cut_off_column_top ! ----------------------------------------------------------------------------- @@ -2486,6 +2517,7 @@ subroutine MOM_state_init_tests(G, GV, tv) real, dimension(nk+1) :: e integer :: k real :: P_tot, P_t, P_b, z_out + type(remapping_CS), pointer :: remap_CS => NULL() do k = 1, nk h(k) = 100. @@ -2521,7 +2553,7 @@ subroutine MOM_state_init_tests(G, GV, tv) write(0,*) '' write(0,*) h call cut_off_column_top(nk, tv, GV%Rho0, G%G_earth, -e(nk+1), GV%Angstrom, & - T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h) + T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS) write(0,*) h end subroutine MOM_state_init_tests diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 1a08d800f2..f4f8e4d326 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -79,7 +79,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) real :: wc ! half-width of the trough real :: ly ! domain width (across ice flow) real :: bx, by, xtil ! dummy vatiables - + logical :: is_2D ! If true, use 2D setup + ! G%ieg and G%jeg are the last indices in the global domain ! real, dimension (G%ieg) :: xtil ! eq. 3 ! real, dimension (G%ieg) :: bx ! eq. 2 @@ -99,23 +100,44 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) call log_version(param_file, mod, version, "") call get_param(param_file, mod, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0) + call get_param(param_file, mod, "ISOMIP_2D",is_2D,'If true, use a 2D setup.', default=.false.) ! The following variables should be transformed into runtime parameters? bmax=720.0; b0=-150.0; b2=-728.8; b4=343.91; b6=-50.57 xbar=300.0E3; dc=500.0; fc=4.0E3; wc=24.0E3; ly=80.0E3 bx = 0.0; by = 0.0; xtil = 0.0 - do j=js,je ; do i=is,ie - xtil = G%geoLonT(i,j)*1.0e3/xbar - bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & - (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) -! depth is positive - D(i,j) = -max(bx+by,-bmax) - if (D(i,j) > max_depth) D(i,j) = max_depth - if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth - enddo ; enddo + if (is_2D) then + do j=js,je ; do i=is,ie + ! 2D setup + xtil = G%geoLonT(i,j)*1.0e3/xbar + !xtil = 450*1.0e3/xbar + bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 + !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & + ! (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + + by = (dc/(1.+exp(-2.*(20.*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(20.*1.0e3- ly/2. + wc)/fc))) + + D(i,j) = -max(bx+by,-bmax) + if (D(i,j) > max_depth) D(i,j) = max_depth + if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth + enddo ; enddo + + else + do j=js,je ; do i=is,ie + ! 3D + xtil = G%geoLonT(i,j)*1.0e3/xbar + bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6 + by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + & + (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc))) + + D(i,j) = -max(bx+by,-bmax) + if (D(i,j) > max_depth) D(i,j) = max_depth + if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth + enddo ; enddo + endif end subroutine ISOMIP_initialize_topography ! ----------------------------------------------------------------------------- @@ -132,7 +154,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers !! to any available thermodynamic !! fields, including eq. of state. - + ! Local variables real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! @@ -154,18 +176,18 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Initial thicknesses for isopycnal coordinates - call get_param(param_file, mod, "T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9) - call get_param(param_file, mod, "S_SUF", s_sur, 'Salinity at the surface (interface)', default=33.8) - call get_param(param_file, mod, "T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9) - call get_param(param_file, mod, "S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55) + call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9) + call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8) + call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9) + call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55) ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) - write (*,*)'Surface density is:', rho_sur + !write (*,*)'Surface density is:', rho_sur call calculate_density(t_bot,s_bot,0.0,rho_bot,tv%eqn_of_state) - write (*,*)'Bottom density is:', rho_bot + !write (*,*)'Bottom density is:', rho_bot rho_range = rho_bot - rho_sur - write (*,*)'Density range is:', rho_range + !write (*,*)'Density range is:', rho_range ! Construct notional interface positions e0(1) = 0. @@ -173,6 +195,8 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range e0(k) = min( 0., e0(k) ) ! Bound by surface e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model + !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) + enddo e0(nz+1) = -G%max_depth @@ -233,39 +257,34 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, integer :: i, j, k, is, ie, js, je, nz real :: x, ds, dt, rho_sur, rho_bot real :: xi0, xi1, dxi, r, S_sur, T_sur, S_bot, T_bot, S_range, T_range -! real :: T_ref, T_Light, T_Dense, S_ref, S_Light, S_Dense, a1, frac_dense, k_frac, res_rat, k_light real :: z ! vertical position in z space character(len=40) :: verticalCoordinate, density_profile -real :: rho_light, delta_rho, z_tmp, rho_tmp + real :: rho_tmp is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE) - call get_param(param_file, mod, "T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9,& + call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9,& do_not_log=.true.) - call get_param(param_file, mod, "S_SUF", s_sur, 'Salinity at the surface (interface)', default=33.8,& + call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8,& do_not_log=.true.) - call get_param(param_file, mod, "T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9,& + call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9,& do_not_log=.true.) - call get_param(param_file, mod, "S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55,& + call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55,& do_not_log=.true.) call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state) - write (*,*)'Density in the surface layer:', rho_sur + !write (*,*)'Density in the surface layer:', rho_sur call calculate_density(t_bot,s_bot,0.0,rho_bot,eqn_of_state) - write (*,*)'Density in the bottom layer::', rho_bot - - call get_param(param_file, mod, "LIGHTEST_DENSITY", rho_light) - call get_param(param_file, mod, "DENSITY_RANGE", delta_rho) + !write (*,*)'Density in the bottom layer::', rho_bot S_range = s_sur - s_bot T_range = t_sur - t_bot - write(*,*)'s_bot, s_sur, S_range, t_bot, t_sur, T_range', s_bot, s_sur, S_range, t_bot, t_sur, T_range select case ( coordinateMode(verticalCoordinate) ) - case ( REGRIDDING_SIGMA, REGRIDDING_ZSTAR, REGRIDDING_RHO ) + case ( REGRIDDING_RHO, REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) S_range = S_range / G%max_depth ! Convert S_range into dS/dz T_range = T_range / G%max_depth ! Convert T_range into dT/dz do j=js,je ; do i=is,ie @@ -277,11 +296,24 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer enddo enddo ; enddo -i=G%iec; j=G%jec -do k = 1,nz - call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) - write(0,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) -enddo + +!i=G%iec; j=G%jec +!do k = 1,nz +! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state) +! !write(*,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) +!enddo + +! case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) +! +! do j=js,je ; do i=is,ie +! xi0 = 0.0; +! do k = 1,nz +! xi1 = xi0 + h(i,j,k) / G%max_depth; +! S(i,j,k) = S_sur + 0.5 * S_range * (xi0 + xi1); +! T(i,j,k) = T_sur + 0.5 * T_range * (xi0 + xi1); +! xi0 = xi1; +! enddo +! enddo ; enddo case default call MOM_error(FATAL,"isomip_initialize: "// & @@ -307,57 +339,123 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) !! parameter values. type(ALE_sponge_CS), pointer :: CSp !< A pointer that is set to point to !! the control structure for this module. - -! real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! A temporary array for eta. + ! Local variables real :: T(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for temp real :: S(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for salt real :: RHO(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for RHO real :: h(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for thickness real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate, in s-1. - real :: S_ref, T_ref; ! Reference salinity and temerature within - ! surface layer - real :: S_range, T_range; ! Range of salinities and temperatures over the - ! vertical - + real :: TNUDG ! Nudging time scale, days + real :: S_sur, T_sur; ! Surface salinity and temerature in sponge + real :: S_bot, T_bot; ! Bottom salinity and temerature in sponge + real :: t_ref, s_ref ! reference T and S + real :: rho_sur, rho_bot, rho_range, t_range, s_range real :: e0(SZK_(G)) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! real :: eta1D(SZK_(G)+1) ! Interface height relative to the sea surface ! ! positive upward, in m. - real :: min_depth, dummy1, z, Bmax - real :: damp, rho_dummy + real :: min_depth, dummy1, z, delta_h + real :: damp, rho_dummy, min_thickness, rho_tmp, xi0 + character(len=40) :: verticalCoordinate + character(len=40) :: mod = "ISOMIP_initialize_sponges" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + call get_param(PF,mod,"MIN_THICKNESS",min_thickness,'Minimum layer thickness',units='m',default=1.e-3) + + call get_param(PF,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & + default=DEFAULT_COORDINATE_MODE) + + call get_param(PF, mod, "ISOMIP_TNUDG", TNUDG, 'Nudging time scale for sponge layers (days)', default=0.0) + + call get_param(PF, mod, "T_REF", t_ref, 'Reference temperature', default=10.0,& + do_not_log=.true.) + + call get_param(PF, mod, "S_REF", s_ref, 'Reference salinity', default=35.0,& + do_not_log=.true.) + + call get_param(PF, mod, "ISOMIP_S_SUR_SPONGE", s_sur, 'Surface salinity in sponge layer.', default=s_ref) + + call get_param(PF, mod, "ISOMIP_S_BOT_SPONGE", s_bot, 'Bottom salinity in sponge layer.', default=s_ref) + + call get_param(PF, mod, "ISOMIP_T_SUR_SPONGE", t_sur, 'Surface temperature in sponge layer.', default=t_ref) + + call get_param(PF, mod, "ISOMIP_T_BOT_SPONGE", t_bot, 'Bottom temperature in sponge layer.', default=t_ref) + T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0; RHO(:,:,:) = 0.0 - S_ref = 33.8; S_range = 0.90 - T_ref = -1.9; T_range = 2.9 - Bmax = 720.0 + S_range = s_sur - s_bot + T_range = t_sur - t_bot ! Set up sponges for ISOMIP configuration call get_param(PF, mod, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0) -! GMM: set thickness, I am using same config. as the initial thickness for now - do k=1,nz - e0(k) = -G%max_depth * real(k-1) / real(nz) - enddo - - do j=js,je ; do i=is,ie - eta1D(nz+1) = -1.0*G%bathyT(i,j) - do k=nz,1,-1 - eta1D(k) = e0(k) - if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then - eta1D(k) = eta1D(k+1) + GV%Angstrom_z - h(i,j,k) = GV%Angstrom_z - else - h(i,j,k) = eta1D(k) - eta1D(k+1) - endif - enddo - enddo; enddo + + select case ( coordinateMode(verticalCoordinate) ) + + case ( REGRIDDING_LAYER, REGRIDDING_RHO ) + ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT + call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) + !write (*,*)'Surface density in sponge:', rho_sur + call calculate_density(t_bot,s_bot,0.0,rho_bot,tv%eqn_of_state) + !write (*,*)'Bottom density in sponge:', rho_bot + rho_range = rho_bot - rho_sur + !write (*,*)'Density range in sponge:', rho_range + + ! Construct notional interface positions + e0(1) = 0. + do K=2,nz + e0(k) = -G%max_depth * ( 0.5 * ( GV%Rlay(k-1) + GV%Rlay(k) ) - rho_sur ) / rho_range + e0(k) = min( 0., e0(k) ) ! Bound by surface + e0(k) = max( -G%max_depth, e0(k) ) ! Bound by possible deepest point in model + !write(*,*)'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k) + + enddo + e0(nz+1) = -G%max_depth + + ! Calculate thicknesses + do j=js,je ; do i=is,ie + eta1D(nz+1) = -1.0*G%bathyT(i,j) + do k=nz,1,-1 + eta1D(k) = e0(k) + if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_z)) then + eta1D(k) = eta1D(k+1) + GV%Angstrom_z + h(i,j,k) = GV%Angstrom_z + else + h(i,j,k) = eta1D(k) - eta1D(k+1) + endif + enddo + enddo ; enddo + + case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates + do j=js,je ; do i=is,ie + eta1D(nz+1) = -1.0*G%bathyT(i,j) + do k=nz,1,-1 + eta1D(k) = -G%max_depth * real(k-1) / real(nz) + if (eta1D(k) < (eta1D(k+1) + min_thickness)) then + eta1D(k) = eta1D(k+1) + min_thickness + h(i,j,k) = min_thickness + else + h(i,j,k) = eta1D(k) - eta1D(k+1) + endif + enddo + enddo ; enddo + + case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates + do j=js,je ; do i=is,ie + delta_h = G%bathyT(i,j) / dfloat(nz) + h(i,j,:) = delta_h + end do ; end do + + case default + call MOM_error(FATAL,"ISOMIP_initialize_sponges: "// & + "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") + + end select ! Here the inverse damping time, in s-1, is set. Set Idamp to 0 ! ! wherever there is no sponge, and the subroutines that are called ! @@ -369,8 +467,7 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) ! 1 / day dummy1=(G%geoLonT(i,j)-790.0)/(800.0-790.0) -! damp = 1.0/10.0 * max(0.0,dummy1) - damp = 1.0/0.04166 * max(0.0,dummy1) ! one hour + damp = 1.0/TNUDG * max(0.0,dummy1) else ; damp=0.0 endif @@ -395,24 +492,37 @@ subroutine ISOMIP_initialize_sponges(G,GV, tv, PF, CSp) call initialize_ALE_sponge(Idamp,h, nz, G, PF, CSp) ! setup temp and salt at the sponge zone + select case ( coordinateMode(verticalCoordinate) ) + + case ( REGRIDDING_SIGMA, REGRIDDING_ZSTAR, REGRIDDING_RHO ) + S_range = S_range / G%max_depth ! Convert S_range into dS/dz + T_range = T_range / G%max_depth ! Convert T_range into dT/dz do j=js,je ; do i=is,ie + xi0 = -G%bathyT(i,j); do k = nz,1,-1 -! z = (G%bathyT(i,j)/(nz-1))* (k -1) - z = (G%max_depth/(nz-1))* (k -1) - S(i,j,k) = S_REF + (S_RANGE*z/Bmax) - T(i,j,k) = T_REF + (T_RANGE*z/Bmax) -! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_dummy,tv%eqn_of_state) -! RHO(i,j,k) = rho_dummy + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth in middle of layer + S(i,j,k) = S_sur + S_range * xi0 + T(i,j,k) = T_sur + T_range * xi0 + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer enddo enddo ; enddo + i=G%iec; j=G%jec + do k = 1,nz + call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state) + !write(*,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k) + enddo + + case default + call MOM_error(FATAL,"isomip_initialize_sponge: "// & + "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE") + + end select + ! Now register all of the fields which are damped in the sponge. ! ! By default, momentum is advected vertically within the sponge, but ! ! momentum is typically not damped within the sponge. ! -! At this point, the ISOMIP configuration is done. The following are here as a -! template for other configurations. - ! The remaining calls to set_up_sponge_field can be in any order. ! if ( associated(tv%T) ) then call set_up_ALE_sponge_field(T,G,tv%T,CSp)