Skip to content

Commit

Permalink
Merge pull request #53 from NOAA-GFDL/dev/gfdl
Browse files Browse the repository at this point in the history 
Merge in dev/gfdl updates
  • Loading branch information
@wrongkindofdoctor
wrongkindofdoctor authored Mar 30, 2020
2 parents c245357 + 02f257b commit 86c2a7d
Show file tree
Hide file tree
Showing 43 changed files with 552 additions and 524 deletions.
70 changes: 35 additions & 35 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ module MOM_regridding
!> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2].
real :: min_thickness

!> Reference pressure for potential density calculations (Pa)
!> Reference pressure for potential density calculations [Pa]
real :: ref_pressure = 2.e7

!> Weight given to old coordinate when blending between new and old grids [nondim]
Expand Down Expand Up @@ -203,6 +203,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
real :: dz_fixed_sfc, Rho_avg_depth, nlay_sfc_int
real :: adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha
real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3]
integer :: nz_fixed_sfc, k, nzf(4)
real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m]
! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
Expand Down Expand Up @@ -559,7 +560,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
call get_param(param_file, mdl, "HALOCLINE_FILTER_LENGTH", filt_len, &
"A length scale over which to smooth the temperature and "//&
"salinity before identifying erroneously unstable haloclines.", &
units="m", default=2.0)
units="m", default=2.0, scale=GV%m_to_H)
call get_param(param_file, mdl, "HALOCLINE_STRAT_TOL", strat_tol, &
"A tolerance for the ratio of the stratification of the "//&
"apparent coordinate stratification to the actual value "//&
Expand All @@ -574,26 +575,26 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m

if (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adaptTimeRatio, &
"Ratio of ALE timestep to grid timescale.", units="nondim", default=1e-1)
"Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1)
call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptZoom, &
"Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H)
"Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H)
call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptZoomCoeff, &
"Coefficient of near-surface zooming diffusivity.", &
units="nondim", default=0.2)
"Coefficient of near-surface zooming diffusivity.", units="nondim", default=0.2)
call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptBuoyCoeff, &
"Coefficient of buoyancy diffusivity.", &
units="nondim", default=0.8)
"Coefficient of buoyancy diffusivity.", units="nondim", default=0.8)
call get_param(param_file, mdl, "ADAPT_ALPHA", adaptAlpha, &
"Scaling on optimization tendency.", &
units="nondim", default=1.0)
"Scaling on optimization tendency.", units="nondim", default=1.0)
call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmpLogical, &
"If true, make a HyCOM-like mixed layer by preventing interfaces "//&
"from being shallower than the depths specified by the regridding coordinate.", &
default=.false.)
"If true, make a HyCOM-like mixed layer by preventing interfaces "//&
"from being shallower than the depths specified by the regridding coordinate.", &
default=.false.)
call get_param(param_file, mdl, "ADAPT_DRHO0", adaptDrho0, &
"Reference density difference for stratification-dependent diffusion.", &
units="kg m-3", default=0.5, scale=US%kg_m3_to_R)

call set_regrid_params(CS, adaptTimeRatio=adaptTimeRatio, adaptZoom=adaptZoom, &
adaptZoomCoeff=adaptZoomCoeff, adaptBuoyCoeff=adaptBuoyCoeff, adaptAlpha=adaptAlpha, &
adaptDoMin=tmpLogical)
adaptDoMin=tmpLogical, adaptDrho0=US%R_to_kg_m3*adaptDrho0)
endif

if (main_parameters .and. coord_is_state_dependent) then
Expand Down Expand Up @@ -1012,9 +1013,9 @@ end subroutine check_grid_column
subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
type(regridding_CS), intent(in) :: CS !< Regridding control structure
integer, intent(in) :: nk !< Number of cells in source grid
real, dimension(nk+1), intent(in) :: z_old !< Old grid position [m]
real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position [m]
real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [m]
real, dimension(nk+1), intent(in) :: z_old !< Old grid position [H ~> m or kg m-2]
real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position [H ~> m or kg m-2]
real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions [H ~> m or kg m-2]
! Local variables
real :: sgn ! The sign convention for downward.
real :: dz_tgt, zr1, z_old_k
Expand Down Expand Up @@ -1156,26 +1157,21 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
!! [H ~> m or kg m-2].
real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage.
real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage [nondim].
! Local variables
integer :: i, j, k
integer :: nz
real :: nominalDepth, totalThickness, dh
real, dimension(SZK_(GV)+1) :: zOld, zNew
real :: minThickness
real :: nominalDepth, totalThickness, dh ! Depths and thicknesses [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: zOld, zNew ! Coordinate interface heights [H ~> m or kg m-2]
integer :: i, j, k, nz
logical :: ice_shelf

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

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

Expand Down Expand Up @@ -1218,7 +1214,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)
#ifdef __DO_SAFETY_CHECKS__
dh=max(nominalDepth,totalThickness)
if (abs(zNew(1)-zOld(1))>(nz-1)*0.5*epsilon(dh)*dh) then
write(0,*) 'min_thickness=',minThickness
write(0,*) 'min_thickness=',CS%min_thickness
write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness
write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz
do k=1,nz+1
Expand Down Expand Up @@ -1350,10 +1346,11 @@ subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
! Local variables
integer :: nz
integer :: i, j, k
real :: nominalDepth, totalThickness
real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: zOld, zNew
real :: h_neglect, h_neglect_edge
real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
#ifdef __DO_SAFETY_CHECKS__
real :: totalThickness
real :: dh
#endif

Expand Down Expand Up @@ -1539,8 +1536,8 @@ subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt
! current interface positions and after tendency term is applied
! positive downward
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt
real, dimension(SZK_(GV)+1) :: zNext
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt ! Interface depths [H ~> m or kg m-2]
real, dimension(SZK_(GV)+1) :: zNext ! New interface depths [H ~> m or kg m-2]

nz = GV%ke

Expand Down Expand Up @@ -2231,7 +2228,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, &
nlay_ML_to_interior, fix_haloclines, halocline_filt_len, &
halocline_strat_tol, integrate_downward_for_e, remap_answers_2018, &
adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin, adaptDrho0)
type(regridding_CS), intent(inout) :: CS !< Regridding control structure
logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the
Expand All @@ -2250,7 +2247,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
!! resolved stratification [nondim]
logical, optional, intent(in) :: fix_haloclines !< Detect regions with much weaker stratification in the coordinate
real, optional, intent(in) :: halocline_filt_len !< Length scale over which to filter T & S when looking for
!! spuriously unstable water mass profiles [m]
!! spuriously unstable water mass profiles [H ~> m or kg m-2]
real, optional, intent(in) :: halocline_strat_tol !< Value of the stratification ratio that defines a problematic
!! halocline region.
logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward
Expand All @@ -2266,6 +2263,8 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
logical, optional, intent(in) :: adaptDoMin !< If true, make a HyCOM-like mixed layer by
!! preventing interfaces from being shallower than
!! the depths specified by the regridding coordinate.
real, optional, intent(in) :: adaptDrho0 !< Reference density difference for stratification-dependent
!! diffusion. [kg m-3]

if (present(interp_scheme)) call set_interp_scheme(CS%interp_CS, interp_scheme)
if (present(boundary_extrapolation)) call set_interp_extrap(CS%interp_CS, boundary_extrapolation)
Expand Down Expand Up @@ -2322,6 +2321,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
if (present(adaptBuoyCoeff)) call set_adapt_params(CS%adapt_CS, adaptBuoyCoeff=adaptBuoyCoeff)
if (present(adaptAlpha)) call set_adapt_params(CS%adapt_CS, adaptAlpha=adaptAlpha)
if (present(adaptDoMin)) call set_adapt_params(CS%adapt_CS, adaptDoMin=adaptDoMin)
if (present(adaptDrho0)) call set_adapt_params(CS%adapt_CS, adaptDrho0=adaptDrho0)
end select

end subroutine set_regrid_params
Expand Down
46 changes: 22 additions & 24 deletions src/ALE/coord_hycom.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ module coord_hycom
!> Number of layers/levels in generated grid
integer :: nk

!> Nominal near-surface resolution
!> Nominal near-surface resolution [Z ~> m]
real, allocatable, dimension(:) :: coordinateResolution

!> Nominal density of interfaces [R ~> kg m-3]
Expand All @@ -24,10 +24,10 @@ module coord_hycom
!> Density scaling factor [R m3 kg-1 ~> 1]
real :: kg_m3_to_R

!> Maximum depths of interfaces
!> Maximum depths of interfaces [H ~> m or kg m-2]
real, allocatable, dimension(:) :: max_interface_depths

!> Maximum thicknesses of layers
!> Maximum thicknesses of layers [H ~> m or kg m-2]
real, allocatable, dimension(:) :: max_layer_thickness

!> Interpolation control structure
Expand All @@ -42,7 +42,7 @@ module coord_hycom
subroutine init_coord_hycom(CS, nk, coordinateResolution, target_density, interp_CS, rho_scale)
type(hycom_CS), pointer :: CS !< Unassociated pointer to hold the control structure
integer, intent(in) :: nk !< Number of layers in generated grid
real, dimension(nk), intent(in) :: coordinateResolution !< Nominal near-surface resolution [m]
real, dimension(nk), intent(in) :: coordinateResolution !< Nominal near-surface resolution [Z ~> m]
real, dimension(nk+1),intent(in) :: target_density !< Interface target densities [R ~> kg m-3]
type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation
real, optional, intent(in) :: rho_scale !< A dimensional scaling factor for target_density
Expand Down Expand Up @@ -76,8 +76,8 @@ end subroutine end_coord_hycom
!> This subroutine can be used to set the parameters for the coord_hycom module
subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, interp_CS)
type(hycom_CS), pointer :: CS !< Coordinate control structure
real, dimension(:), optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces in m
real, dimension(:), optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers in m
real, dimension(:), optional, intent(in) :: max_interface_depths !< Maximum depths of interfaces [H ~> m or kg m-2]
real, dimension(:), optional, intent(in) :: max_layer_thickness !< Maximum thicknesses of layers [H ~> m or kg m-2]
type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation

if (.not. associated(CS)) call MOM_error(FATAL, "set_hycom_params: CS not associated")
Expand All @@ -102,33 +102,31 @@ end subroutine set_hycom_params
!> Build a HyCOM coordinate column
subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, &
z_col, z_col_new, zScale, h_neglect, h_neglect_edge)
type(hycom_CS), intent(in) :: CS !< Coordinate control structure
type(hycom_CS), intent(in) :: CS !< Coordinate control structure
type(EOS_type), pointer :: eqn_of_state !< Equation of state structure
integer, intent(in) :: nz !< Number of levels
integer, intent(in) :: nz !< Number of levels
real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2])
real, dimension(nz), intent(in) :: T !< Temperature of column [degC]
real, dimension(nz), intent(in) :: S !< Salinity of column [ppt]
real, dimension(nz), intent(in) :: h !< Layer thicknesses, in [m] or [H ~> m or kg m-2]
real, dimension(nz), intent(in) :: T !< Temperature of column [degC]
real, dimension(nz), intent(in) :: S !< Salinity of column [ppt]
real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, dimension(nz), intent(in) :: p_col !< Layer pressure [Pa]
real, dimension(nz+1), intent(in) :: z_col !< Interface positions relative to the surface [H ~> m or kg m-2]
real, dimension(CS%nk+1), intent(inout) :: z_col_new !< Absolute positions of interfaces
real, optional, intent(in) :: zScale !< Scaling factor from the input thicknesses in [m]
!! to desired units for zInterface, perhaps m_to_H.
real, optional, intent(in) :: h_neglect !< A negligibly small width for the
!! purpose of cell reconstructions
!! in the same units as h.
real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
!! for the purpose of edge value calculations
!! in the same units as h0.
real, dimension(CS%nk+1), intent(inout) :: z_col_new !< Absolute positions of interfaces [H ~> m or kg m-2]
real, optional, intent(in) :: zScale !< Scaling factor from the input coordinate thicknesses in [Z ~> m]
!! to desired units for zInterface, perhaps GV%Z_to_H.
real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose of
!! cell reconstruction [H ~> m or kg m-2]
real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose of
!! edge value calculation [H ~> m or kg m-2]

! Local variables
integer :: k
real, dimension(nz) :: rho_col ! Layer densities in a column [R ~> kg m-3]
real, dimension(CS%nk) :: h_col_new ! New layer thicknesses
real :: z_scale
real :: stretching ! z* stretching, converts z* to z.
real :: nominal_z ! Nominal depth of interface when using z* [Z ~> m]
real :: hNew
real :: z_scale ! A scaling factor from the input thicknesses to the target thicknesses,
! perhaps 1 or a factor in [H Z-1 ~> 1 or kg m-3]
real :: stretching ! z* stretching, converts z* to z [nondim].
real :: nominal_z ! Nominal depth of interface when using z* [H ~> m or kg m-2]
logical :: maximum_depths_set ! If true, the maximum depths of interface have been set.
logical :: maximum_h_set ! If true, the maximum layer thicknesses have been set.

Expand Down
Loading

0 comments on commit 86c2a7d

Please sign in to comment.