Skip to content

Commit

Permalink
Merge pull request #56 from NOAA-GFDL/dev/gfdl
Browse files Browse the repository at this point in the history
Merge in latest dev/gfdl updates
  • Loading branch information
wrongkindofdoctor authored May 1, 2020
2 parents e9281f4 + 6603451 commit 496617c
Show file tree
Hide file tree
Showing 78 changed files with 3,876 additions and 3,702 deletions.
61 changes: 30 additions & 31 deletions config_src/coupled_driver/MOM_surface_forcing_gfdl.F90
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,8 @@ module MOM_surface_forcing_gfdl
real :: latent_heat_fusion !< Latent heat of fusion [J kg-1]
real :: latent_heat_vapor !< Latent heat of vaporization [J kg-1]

real :: max_p_surf !< The maximum surface pressure that can be
!! exerted by the atmosphere and floating sea-ice [Pa].
real :: max_p_surf !< The maximum surface pressure that can be exerted by
!! the atmosphere and floating sea-ice [R L2 T-2 ~> Pa].
!! This is needed because the FMS coupling structure
!! does not limit the water that can be frozen out
!! of the ocean and the ice-ocean heat fluxes are
Expand Down Expand Up @@ -101,12 +101,12 @@ module MOM_surface_forcing_gfdl

logical :: rigid_sea_ice !< If true, sea-ice exerts a rigidity that acts to damp surface
!! deflections (especially surface gravity waves). The default is false.
real :: G_Earth !< Gravitational acceleration [m s-2]
real :: Kv_sea_ice !< Viscosity in sea-ice that resists sheared vertical motions [m2 s-1]
real :: density_sea_ice !< Typical density of sea-ice [kg m-3]. The value is only used to convert
real :: g_Earth !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
real :: Kv_sea_ice !< Viscosity in sea-ice that resists sheared vertical motions [L4 Z-2 T-1 ~> m2 s-1]
real :: density_sea_ice !< Typical density of sea-ice [R ~> kg m-3]. The value is only used to convert
!! the ice pressure into appropriate units for use with Kv_sea_ice.
real :: rigid_sea_ice_mass !< A mass per unit area of sea-ice beyond which sea-ice viscosity
!! becomes effective [kg m-2], typically of order 1000 kg m-2.
!! becomes effective [R Z ~> kg m-2], typically of order 1000 kg m-2.
logical :: allow_flux_adjustments !< If true, use data_override to obtain flux adjustments

logical :: restore_salt !< If true, the coupled MOM driver adds a term to restore surface
Expand Down Expand Up @@ -466,7 +466,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
endif

if (associated(IOB%mass_berg)) then
fluxes%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j)
fluxes%mass_berg(i,j) = US%m_to_Z*US%kg_m3_to_R * IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%mass_berg(i-i0,j-j0), G%mask2dT(i,j), i, j, 'mass_berg', G)
endif
Expand Down Expand Up @@ -548,14 +548,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
if (associated(IOB%p)) then
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%p(i-i0,j-j0), G%mask2dT(i,j), i, j, 'p', G)
enddo ; enddo
else
do j=js,je ; do i=is,ie
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%p(i-i0,j-j0), G%mask2dT(i,j), i, j, 'p', G)
Expand Down Expand Up @@ -669,14 +669,14 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: &
rigidity_at_h, & ! Ice rigidity at tracer points [m3 s-1]
rigidity_at_h, & ! Ice rigidity at tracer points [L4 Z-1 T-1 ~> m3 s-1]
net_mass_src, & ! A temporary of net mass sources [kg m-2 s-1].
ustar_tmp ! A temporary array of ustar values [Z T-1 ~> m s-1].

real :: I_GEarth ! 1.0 / G_Earth [s2 m-1]
real :: Kv_rho_ice ! (CS%kv_sea_ice / CS%density_sea_ice) [m5 s-1 kg-1]
real :: mass_ice ! mass of sea ice at a face [kg m-2]
real :: mass_eff ! effective mass of sea ice for rigidity [kg m-2]
real :: I_GEarth ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-1]
real :: Kv_rho_ice ! (CS%Kv_sea_ice / CS%density_sea_ice) [L4 Z-2 T-1 R-1 ~> m5 s-1 kg-1]
real :: mass_ice ! mass of sea ice at a face [R Z ~> kg m-2]
real :: mass_eff ! effective mass of sea ice for rigidity [R Z ~> kg m-2]
real :: wt1, wt2 ! Relative weights of previous and current values of ustar, ND.

integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
Expand Down Expand Up @@ -751,12 +751,12 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
if (associated(IOB%p)) then
if (CS%max_p_surf >= 0.0) then
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
forces%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = MIN(forces%p_surf_full(i,j),CS%max_p_surf)
enddo ; enddo
else
do j=js,je ; do i=is,ie
forces%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0)
forces%p_surf_full(i,j) = G%mask2dT(i,j) * US%kg_m3_to_R*US%m_s_to_L_T**2*IOB%p(i-i0,j-j0)
forces%p_surf(i,j) = forces%p_surf_full(i,j)
enddo ; enddo
endif
Expand Down Expand Up @@ -816,13 +816,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
enddo ; enddo ; endif

if (associated(IOB%mass_berg)) then ; do j=js,je ; do i=is,ie
forces%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j)
forces%mass_berg(i,j) = US%m_to_Z*US%kg_m3_to_R * IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j)
enddo ; enddo ; endif

! Obtain sea ice related dynamic fields
if (associated(IOB%ice_rigidity)) then
do j=js,je ; do i=is,ie
rigidity_at_h(i,j) = IOB%ice_rigidity(i-i0,j-j0) * G%mask2dT(i,j)
rigidity_at_h(i,j) = US%m_to_L**3*US%Z_to_L*US%T_to_s * IOB%ice_rigidity(i-i0,j-j0) * G%mask2dT(i,j)
enddo ; enddo
call pass_var(rigidity_at_h, G%Domain, halo=1)
do I=is-1,ie ; do j=js,je
Expand All @@ -837,23 +837,21 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_

if (CS%rigid_sea_ice) then
call pass_var(forces%p_surf_full, G%Domain, halo=1)
I_GEarth = 1.0 / CS%G_Earth
Kv_rho_ice = (CS%kv_sea_ice / CS%density_sea_ice)
I_GEarth = 1.0 / CS%g_Earth
Kv_rho_ice = (CS%Kv_sea_ice / CS%density_sea_ice)
do I=is-1,ie ; do j=js,je
mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i+1,j)) * I_GEarth
mass_eff = 0.0
if (mass_ice > CS%rigid_sea_ice_mass) then
mass_eff = (mass_ice - CS%rigid_sea_ice_mass) **2 / &
(mass_ice + CS%rigid_sea_ice_mass)
mass_eff = (mass_ice - CS%rigid_sea_ice_mass)**2 / (mass_ice + CS%rigid_sea_ice_mass)
endif
forces%rigidity_ice_u(I,j) = forces%rigidity_ice_u(I,j) + Kv_rho_ice * mass_eff
enddo ; enddo
do i=is,ie ; do J=js-1,je
mass_ice = min(forces%p_surf_full(i,j), forces%p_surf_full(i,j+1)) * I_GEarth
mass_eff = 0.0
if (mass_ice > CS%rigid_sea_ice_mass) then
mass_eff = (mass_ice - CS%rigid_sea_ice_mass) **2 / &
(mass_ice + CS%rigid_sea_ice_mass)
mass_eff = (mass_ice - CS%rigid_sea_ice_mass)**2 / (mass_ice + CS%rigid_sea_ice_mass)
endif
forces%rigidity_ice_v(i,J) = forces%rigidity_ice_v(i,J) + Kv_rho_ice * mass_eff
enddo ; enddo
Expand Down Expand Up @@ -1299,8 +1297,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS)
"needed because the FMS coupling structure does not "//&
"limit the water that can be frozen out of the ocean and "//&
"the ice-ocean heat fluxes are treated explicitly. No "//&
"limit is applied if a negative value is used.", units="Pa", &
default=-1.0)
"limit is applied if a negative value is used.", &
units="Pa", default=-1.0, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
call get_param(param_file, mdl, "RESTORE_SALINITY", CS%restore_salt, &
"If true, the coupled driver will add a globally-balanced "//&
"fresh-water flux that drives sea-surface salinity "//&
Expand Down Expand Up @@ -1515,18 +1513,19 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS)
if (CS%rigid_sea_ice) then
call get_param(param_file, mdl, "G_EARTH", CS%g_Earth, &
"The gravitational acceleration of the Earth.", &
units="m s-2", default = 9.80)
units="m s-2", default = 9.80, scale=US%Z_to_m*US%m_s_to_L_T**2)
call get_param(param_file, mdl, "SEA_ICE_MEAN_DENSITY", CS%density_sea_ice, &
"A typical density of sea ice, used with the kinematic "//&
"viscosity, when USE_RIGID_SEA_ICE is true.", units="kg m-3", &
default=900.0)
"viscosity, when USE_RIGID_SEA_ICE is true.", &
units="kg m-3", default=900.0, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "SEA_ICE_VISCOSITY", CS%Kv_sea_ice, &
"The kinematic viscosity of sufficiently thick sea ice "//&
"for use in calculating the rigidity of sea ice.", &
units="m2 s-1", default=1.0e9)
units="m2 s-1", default=1.0e9, scale=US%Z_to_L**2*US%m_to_L**2*US%T_to_s)
call get_param(param_file, mdl, "SEA_ICE_RIGID_MASS", CS%rigid_sea_ice_mass, &
"The mass of sea-ice per unit area at which the sea-ice "//&
"starts to exhibit rigidity", units="kg m-2", default=1000.0)
"starts to exhibit rigidity", &
units="kg m-2", default=1000.0, scale=US%kg_m3_to_R*US%m_to_Z)
endif

call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, &
Expand Down
Loading

0 comments on commit 496617c

Please sign in to comment.