Merge branch 'dev/gfdl' into nonBous_CFC_cap

src/ALE/MOM_ALE.F90

@@ -20,7 +20,7 @@ module MOM_ALE
 use MOM_hybgen_unmix,     only : hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix, hybgen_unmix_CS
 use MOM_hybgen_regrid,    only : hybgen_regrid_CS
 use MOM_file_parser,      only : get_param, param_file_type, log_param
-use MOM_interface_heights,only : find_eta
+use MOM_interface_heights,only : find_eta, calc_derived_thermo
 use MOM_open_boundary,    only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
 use MOM_open_boundary,    only : OBC_DIRECTION_N, OBC_DIRECTION_S
 use MOM_regridding,       only : initialize_regridding, regridding_main, end_regridding
@@ -659,6 +659,9 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d
     ! generate new grid
     if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local)
 
+    ! Update the layer specific volumes if necessary
+    if (allocated(tv_local%SpV_avg)) call calc_derived_thermo(tv_local, h, G, GV, US, halo=1)
+
     call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface)
     dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:)
 

src/ALE/MOM_hybgen_unmix.F90

@@ -9,6 +9,7 @@ module MOM_hybgen_unmix
 use MOM_file_parser,     only : get_param, param_file_type, log_param
 use MOM_hybgen_regrid,   only : hybgen_column_init
 use MOM_hybgen_regrid,   only : hybgen_regrid_CS, get_hybgen_regrid_params
+use MOM_interface_heights, only : calc_derived_thermo
 use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chkinv
 use MOM_unit_scaling,    only : unit_scale_type
 use MOM_variables,       only : ocean_grid_type, thermo_var_ptrs
@@ -146,7 +147,8 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
   real :: p_col(GV%ke)      ! A column of reference pressures [R L2 T-2 ~> Pa]
   real :: tracer(GV%ke,max(ntr,1)) ! Columns of each tracer [Conc]
   real :: h_tot             ! Total thickness of the water column [H ~> m or kg m-2]
-  real :: nominalDepth      ! Depth of ocean bottom (positive downward) [H ~> m or kg m-2]
+  real :: dz_tot            ! Vertical distance between the top and bottom of the water column [Z ~> m]
+  real :: nominalDepth      ! Depth of ocean bottom in thickness units (positive downward) [H ~> m or kg m-2]
   real :: h_thin            ! A negligibly small thickness to identify essentially
                             ! vanished layers [H ~> m or kg m-2]
   real :: dilate            ! A factor by which to dilate the target positions from z to z* [nondim]
@@ -169,6 +171,15 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
   h_thin = 1e-6*GV%m_to_H
   debug_conservation = .false. !  Set this to true for debugging
 
+  if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then
+    if (tv%valid_SpV_halo < 0) then
+      mesg = "invalid values of SpV_avg."
+    else
+      mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed."
+    endif
+    call MOM_error(FATAL, "hybgen_unmix called in fully non-Boussinesq mode with "//trim(mesg))
+  endif
+
   p_col(:) = CS%ref_pressure
 
   do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 ; if (G%mask2dT(i,j)>0.) then
@@ -203,13 +214,27 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
     endif
 
     ! The following block of code is used to trigger z* stretching of the targets heights.
-    nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H
-    if (h_tot <= CS%min_dilate*nominalDepth) then
-      dilate = CS%min_dilate
-    elseif (h_tot >= CS%max_dilate*nominalDepth) then
-      dilate = CS%max_dilate
+    if (allocated(tv%SpV_avg)) then  ! This is the fully non-Boussiesq version
+      dz_tot = 0.0
+      do k=1,nk
+        dz_tot = dz_tot + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h_col(k)
+      enddo
+      if (dz_tot <= CS%min_dilate*(G%bathyT(i,j)+G%Z_ref)) then
+        dilate = CS%min_dilate
+      elseif (dz_tot >= CS%max_dilate*(G%bathyT(i,j)+G%Z_ref)) then
+        dilate = CS%max_dilate
+      else
+        dilate = dz_tot / (G%bathyT(i,j)+G%Z_ref)
+      endif
     else
-      dilate = h_tot / nominalDepth
+      nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H
+      if (h_tot <= CS%min_dilate*nominalDepth) then
+        dilate = CS%min_dilate
+      elseif (h_tot >= CS%max_dilate*nominalDepth) then
+        dilate = CS%max_dilate
+      else
+        dilate = h_tot / nominalDepth
+      endif
     endif
 
     terrain_following = (h_tot < dilate*CS%dpns) .and. (CS%dpns >= CS%dsns)
@@ -268,6 +293,9 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
     endif
   endif ; enddo ; enddo !i & j.
 
+  ! Update the layer properties
+  if (allocated(tv%SpV_avg)) call calc_derived_thermo(tv, h, G, GV, US, halo=1)
+
 end subroutine hybgen_unmix
 
 

src/ALE/MOM_regridding.F90

@@ -814,20 +814,47 @@ subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, &
 
   ! Local variables
   real :: nom_depth_H(SZI_(G),SZJ_(G))  !< The nominal ocean depth at each point in thickness units [H ~> m or kg m-2]
+  real :: tot_h(SZI_(G),SZJ_(G))  !< The total thickness of the water column [H ~> m or kg m-2]
+  real :: tot_dz(SZI_(G),SZJ_(G)) !< The total distance between the top and bottom of the water column [Z ~> m]
   real :: Z_to_H  ! A conversion factor used by some routines to convert coordinate
                   ! parameters to depth units [H Z-1 ~> nondim or kg m-3]
   real :: trickGnuCompiler
-  integer :: i, j
+  character(len=128) :: mesg    ! A string for error messages
+  integer :: i, j, k
 
   if (present(PCM_cell)) PCM_cell(:,:,:) = .false.
 
   Z_to_H = US%Z_to_m * GV%m_to_H  ! Often this is equivalent to GV%Z_to_H.
-  do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
-    nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * Z_to_H
-    ! Consider using the following instead:
-    ! nom_depth_H(i,j) = max( (G%bathyT(i,j)+G%Z_ref) * Z_to_H , CS%min_nom_depth )
-    ! if (G%mask2dT(i,j)==0.) nom_depth_H(i,j) = 0.0
-  enddo ; enddo
+
+  if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then
+    if (tv%valid_SpV_halo < 0) then
+      mesg = "invalid values of SpV_avg."
+    else
+      mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed."
+    endif
+    call MOM_error(FATAL, "Regridding_main called in fully non-Boussinesq mode with "//trim(mesg))
+  endif
+
+  if (allocated(tv%SpV_avg)) then  ! This is the fully non-Boussinesq case
+    do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      tot_h(i,j) = 0.0 ; tot_dz(i,j) = 0.0
+    enddo ; enddo
+    do k=1,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      tot_h(i,j) = tot_h(i,j) + h(i,j,k)
+      tot_dz(i,j) = tot_dz(i,j) + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h(i,j,k)
+    enddo ; enddo ; enddo
+    do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      if ((tot_dz(i,j) > 0.0) .and. (G%bathyT(i,j)+G%Z_ref > 0.0)) then
+        nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * (tot_h(i,j) / tot_dz(i,j))
+      else
+        nom_depth_H(i,j) = 0.0
+      endif
+    enddo ; enddo
+  else
+    do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
+      nom_depth_H(i,j) = max((G%bathyT(i,j)+G%Z_ref) * Z_to_H, 0.0)
+    enddo ; enddo
+  endif
 
   select case ( CS%regridding_scheme )
 
@@ -1308,12 +1335,12 @@ subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface )
       ! In sigma coordinates, the bathymetric depth is only used as an arbitrary offset that
       ! cancels out when determining coordinate motion, so referencing the column postions to
       ! the surface is perfectly acceptable, but for preservation of previous answers the
-      ! referencing is done relative to the bottom when in Boussinesq mode.
-      ! if (GV%Boussinesq) then
+      ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
+      if (GV%Boussinesq .or. GV%semi_Boussinesq) then
         nominalDepth = nom_depth_H(i,j)
-      ! else
-      !   nominalDepth = totalThickness
-      ! endif
+      else
+        nominalDepth = totalThickness
+      endif
 
       call build_sigma_column(CS%sigma_CS, nominalDepth, totalThickness, zNew)
 
@@ -1436,12 +1463,12 @@ subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS,
       ! In rho coordinates, the bathymetric depth is only used as an arbitrary offset that
       ! cancels out when determining coordinate motion, so referencing the column postions to
       ! the surface is perfectly acceptable, but for preservation of previous answers the
-      ! referencing is done relative to the bottom when in Boussinesq mode.
-      ! if (GV%Boussinesq) then
+      ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
+      if (GV%Boussinesq .or. GV%semi_Boussinesq) then
         nominalDepth = nom_depth_H(i,j)
-      ! else
-      !   nominalDepth = totalThickness
-      ! endif
+      else
+        nominalDepth = totalThickness
+      endif
 
       ! Determine absolute interface positions
       zOld(nz+1) = - nominalDepth

src/core/MOM.F90

@@ -34,7 +34,7 @@ module MOM
 use MOM_error_handler,        only : MOM_set_verbosity, callTree_showQuery
 use MOM_error_handler,        only : callTree_enter, callTree_leave, callTree_waypoint
 use MOM_file_parser,          only : read_param, get_param, log_version, param_file_type
-use MOM_forcing_type,         only : forcing, mech_forcing
+use MOM_forcing_type,         only : forcing, mech_forcing, find_ustar
 use MOM_forcing_type,         only : MOM_forcing_chksum, MOM_mech_forcing_chksum
 use MOM_get_input,            only : Get_MOM_Input, directories
 use MOM_io,                   only : MOM_io_init, vardesc, var_desc
@@ -91,7 +91,7 @@ module MOM
 use MOM_grid,                  only : set_first_direction, rescale_grid_bathymetry
 use MOM_hor_index,             only : hor_index_type, hor_index_init
 use MOM_hor_index,             only : rotate_hor_index
-use MOM_interface_heights,     only : find_eta, calc_derived_thermo
+use MOM_interface_heights,     only : find_eta, calc_derived_thermo, thickness_to_dz
 use MOM_interface_filter,      only : interface_filter, interface_filter_init, interface_filter_end
 use MOM_interface_filter,      only : interface_filter_CS
 use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end
@@ -544,8 +544,13 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
                         ! the end of a stepping cycle (whatever that may mean).
   logical :: therm_reset ! If true, reset running sums of thermodynamic quantities.
   real :: cycle_time    ! The length of the coupled time-stepping cycle [T ~> s].
+  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
+    U_star      ! The wind friction velocity, calculated using the Boussinesq reference density or
+                ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]
   real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
     ssh         ! sea surface height, which may be based on eta_av [Z ~> m]
+  real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%GV)) :: &
+    dz          ! Vertical distance across layers [Z ~> m]
 
   real, dimension(:,:,:), pointer :: &
     u => NULL(), & ! u : zonal velocity component [L T-1 ~> m s-1]
@@ -672,13 +677,18 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
 
     dt = time_interval / real(n_max)
     dt_therm = dt ; ntstep = 1
+
+    if (CS%UseWaves .and. associated(fluxes%ustar)) &
+      call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass, halo=1)
+    if (CS%UseWaves .and. associated(fluxes%tau_mag)) &
+      call pass_var(fluxes%tau_mag, G%Domain, clock=id_clock_pass, halo=1)
+
     if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf
     CS%tv%p_surf => NULL()
     if (CS%use_p_surf_in_EOS .and. associated(fluxes%p_surf)) then
       CS%tv%p_surf => fluxes%p_surf
       if (allocated(CS%tv%SpV_avg)) call pass_var(fluxes%p_surf, G%Domain, clock=id_clock_pass)
     endif
-    if (CS%UseWaves) call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass)
   endif
 
   if (therm_reset) then
@@ -722,12 +732,16 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
     if (CS%UseWaves) then
       ! Update wave information, which is presently kept static over each call to step_mom
       call enable_averages(time_interval, Time_start + real_to_time(US%T_to_s*time_interval), CS%diag)
-      call Update_Stokes_Drift(G, GV, US, Waves, h, forces%ustar, time_interval, do_dyn)
+      call find_ustar(forces, CS%tv, U_star, G, GV, US, halo=1)
+      call thickness_to_dz(h, CS%tv, dz, G, GV, US, halo_size=1)
+      call Update_Stokes_Drift(G, GV, US, Waves, dz, U_star, time_interval, do_dyn)
       call disable_averaging(CS%diag)
     endif
   else ! not do_dyn.
     if (CS%UseWaves) then ! Diagnostics are not enabled in this call.
-      call Update_Stokes_Drift(G, GV, US, Waves, h, fluxes%ustar, time_interval, do_dyn)
+      call find_ustar(fluxes, CS%tv, U_star, G, GV, US, halo=1)
+      call thickness_to_dz(h, CS%tv, dz, G, GV, US, halo_size=1)
+      call Update_Stokes_Drift(G, GV, US, Waves, dz, U_star, time_interval, do_dyn)
     endif
   endif
 
@@ -2089,7 +2103,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
   call get_param(param_file, "MOM", "CALC_RHO_FOR_SEA_LEVEL", CS%calc_rho_for_sea_lev, &
                  "If true, the in-situ density is used to calculate the "//&
                  "effective sea level that is returned to the coupler. If false, "//&
-                 "the Boussinesq parameter RHO_0 is used.", default=.false.)
+                 "the Boussinesq parameter RHO_0 is used.", default=non_Bous)
   call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
                  "If true, Temperature and salinity are used as state "//&
                  "variables.", default=.true.)
@@ -2878,8 +2892,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
     endif
   endif
 
-  ! Allocate any derived equation of state fields.
-  if (use_temperature .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
+  ! Allocate any derived densities or other equation of state derived fields.
+  if (.not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
     allocate(CS%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0)
     CS%tv%valid_SpV_halo = -1  ! This array does not yet have any valid data.
   endif
@@ -3261,7 +3275,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS)
   ! Write initial conditions
   if (CS%write_IC) then
     allocate(restart_CSp_tmp)
-    restart_CSP_tmp = CS%restart_CS
+    restart_CSp_tmp = CS%restart_CS
     call restart_registry_lock(restart_CSp_tmp, unlocked=.true.)
     allocate(z_interface(SZI_(G),SZJ_(G),SZK_(GV)+1))
     call find_eta(CS%h, CS%tv, G, GV, US, z_interface, dZref=G%Z_ref)

src/core/MOM_dynamics_split_RK2.F90

@@ -481,7 +481,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     Use_Stokes_PGF = associated(Waves)
     if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
     if (Use_Stokes_PGF) then
-      call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
+      call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
+      call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
 
       ! We are adding Stokes_PGF to hydrostatic PGF here.  The diag PFu/PFv
       ! will therefore report the sum total PGF and we avoid other
@@ -748,7 +749,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
       Use_Stokes_PGF = associated(Waves)
       if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
       if (Use_Stokes_PGF) then
-        call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
+        call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
+        call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
         if (.not.Waves%Passive_Stokes_PGF) then
           do k=1,nz
             do j=js,je ; do I=Isq,Ieq