Merge branch 'NOAA-GFDL:dev/gfdl' into dev/cefi

config_src/external/drifters/MOM_particles.F90

@@ -28,7 +28,7 @@ subroutine particles_init(parts, Grid, Time, dt, u, v, h)
 end subroutine particles_init
 
 !> The main driver the steps updates particles
-subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger)
+subroutine particles_run(parts, time, uo, vo, ho, tv, dt_adv, use_uh)
   ! Arguments
   type(particles), pointer :: parts !< Container for all types and memory
   type(time_type), intent(in) :: time !< Model time
@@ -40,8 +40,8 @@ subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger)
                                            !! that are used to advect tracers [H L2 ~> m3 or kg]
   real, dimension(:,:,:), intent(in) :: ho !< Ocean layer thickness [H ~> m or kg m-2]
   type(thermo_var_ptrs),  intent(in) :: tv !< structure containing pointers to available thermodynamic fields
+  real, intent(in) :: dt_adv !< timestep for advecting particles [s]
   logical :: use_uh !< Flag for whether u and v are weighted by thickness
-  integer, optional, intent(in) :: stagger !< Flag for whether velocities are staggered
 
 end subroutine particles_run
 

src/core/MOM.F90

@@ -1283,6 +1283,18 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
 
   endif ! -------------------------------------------------- end SPLIT
 
+  if (CS%use_particles .and. CS%do_dynamics .and. (.not. CS%use_uh_particles)) then
+    if (CS%thickness_diffuse_first) call MOM_error(WARNING,"particles_run: "//&
+      "Thickness_diffuse_first is true and use_uh_particles is false. "//&
+      "This is usually a bad combination.")
+    !Run particles using unweighted velocity
+    call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, &
+                       CS%tv, dt, CS%use_uh_particles)
+    call particles_to_z_space(CS%particles, h)
+  endif
+
+
+
   ! Update the model's current to reflect wind-wave growth
   if (Waves%Stokes_DDT .and. (.not.Waves%Passive_Stokes_DDT)) then
     do J=jsq,jeq ; do i=is,ie
@@ -1368,23 +1380,21 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   endif
   call disable_averaging(CS%diag)
 
+  ! Advance the dynamics time by dt.
+  CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt
+
   if (CS%use_particles .and. CS%do_dynamics .and. CS%use_uh_particles) then
     !Run particles using thickness-weighted velocity
     call particles_run(CS%particles, Time_local, CS%uhtr, CS%vhtr, CS%h, &
-        CS%tv, CS%use_uh_particles)
-  elseif (CS%use_particles .and. CS%do_dynamics) then
-    !Run particles using unweighted velocity
-    call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, &
-                       CS%tv, CS%use_uh_particles)
+        CS%tv, CS%t_dyn_rel_adv, CS%use_uh_particles)
   endif
 
-
-  ! Advance the dynamics time by dt.
-  CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt
   CS%n_dyn_steps_in_adv = CS%n_dyn_steps_in_adv + 1
   if (CS%alternate_first_direction) then
     call set_first_direction(G, MODULO(G%first_direction+1,2))
     CS%first_dir_restart = real(G%first_direction)
+  elseif (CS%use_particles .and. CS%do_dynamics .and. (.not.CS%use_uh_particles)) then
+    call particles_to_k_space(CS%particles, h)
   endif
   CS%t_dyn_rel_thermo = CS%t_dyn_rel_thermo + dt
   if (abs(CS%t_dyn_rel_thermo) < 1e-6*dt) CS%t_dyn_rel_thermo = 0.0

src/core/MOM_density_integrals.F90

@@ -2001,7 +2001,7 @@ end subroutine diagnose_mass_weight_p
 
 !> Find the depth at which the reconstructed pressure matches P_tgt
 subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, &
-                       rho_ref, G_e, EOS, US, P_b, z_out, z_tol)
+                       rho_ref, G_e, EOS, US, P_b, z_out, z_tol, frac_dp_bugfix)
   real,                  intent(in)  :: T_t !< Potential temperature at the cell top [C ~> degC]
   real,                  intent(in)  :: T_b !< Potential temperature at the cell bottom [C ~> degC]
   real,                  intent(in)  :: S_t !< Salinity at the cell top [S ~> ppt]
@@ -2020,6 +2020,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
   real,                  intent(out) :: P_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa]
   real,                  intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m]
   real,                  intent(in)  :: z_tol !< The tolerance in finding z_out [Z ~> m]
+  logical,               intent(in)  :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
 
   ! Local variables
   real :: dp    ! Pressure thickness of the layer [R L2 T-2 ~> Pa]
@@ -2032,7 +2033,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
   GxRho = G_e * rho_ref
 
   ! Anomalous pressure difference across whole cell
-  dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS)
+  dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS, frac_dp_bugfix)
 
   P_b = P_t + dp ! Anomalous pressure at bottom of cell
 
@@ -2063,7 +2064,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
      call MOM_error(FATAL, 'find_depth_of_pressure_in_cell completes too many iterations: '//msg)
     endif
     z_out = z_t + ( z_b - z_t ) * F_guess
-    Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS) - ( P_tgt - P_t )
+    Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS, frac_dp_bugfix) - ( P_tgt - P_t )
 
     if (Pa<Pa_left) then
       write(msg,*) Pa_left,Pa,Pa_right,P_t-P_tgt,P_b-P_tgt
@@ -2119,7 +2120,7 @@ end subroutine avg_specific_vol
 
 !> Returns change in anomalous pressure change from top to non-dimensional
 !! position pos between z_t and z_b [R L2 T-2 ~> Pa]
-real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS)
+real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS, frac_dp_bugfix)
   real,           intent(in)  :: T_t !< Potential temperature at the cell top [C ~> degC]
   real,           intent(in)  :: T_b !< Potential temperature at the cell bottom [C ~> degC]
   real,           intent(in)  :: S_t !< Salinity at the cell top [S ~> ppt]
@@ -2131,6 +2132,7 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO
   real,           intent(in)  :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
   real,           intent(in)  :: pos !< The fractional vertical position, 0 to 1 [nondim]
   type(EOS_type), intent(in)  :: EOS !< Equation of state structure
+  logical,        intent(in)  :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
 
   ! Local variables
   real, parameter :: C1_90 = 1.0/90.0  ! A rational constant [nondim]
@@ -2150,7 +2152,11 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO
     ! Salinity and temperature points are linearly interpolated
     S5(n) = top_weight * S_t + bottom_weight * S_b
     T5(n) = top_weight * T_t + bottom_weight * T_b
-    p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
+    if (frac_dp_bugfix) then
+      p5(n) = (-1) * ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
+    else
+      p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
+    endif !bugfix
   enddo
   call calculate_density(T5, S5, p5, rho5, EOS)
   rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref

src/core/MOM_grid.F90

@@ -190,6 +190,10 @@ module MOM_grid
 
   ! These parameters are run-time parameters that are used during some
   ! initialization routines (but not all)
+  real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related
+                        !! variables like len_lat and len_lon into rescaled horizontal distance
+                        !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or
+                        !! is 0 for a non-Cartesian grid.
   real :: south_lat     !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m]
   real :: west_lon      !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m]
   real :: len_lat       !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m]

src/core/MOM_transcribe_grid.F90

@@ -135,6 +135,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
   ! Copy various scalar variables and strings.
   oG%x_axis_units = dG%x_axis_units ; oG%y_axis_units = dG%y_axis_units
   oG%x_ax_unit_short = dG%x_ax_unit_short ; oG%y_ax_unit_short = dG%y_ax_unit_short
+  oG%grid_unit_to_L = dG%grid_unit_to_L
   oG%areaT_global = dG%areaT_global ; oG%IareaT_global = dG%IareaT_global
   oG%south_lat = dG%south_lat ; oG%west_lon  = dG%west_lon
   oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon
@@ -296,6 +297,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
   ! Copy various scalar variables and strings.
   dG%x_axis_units = oG%x_axis_units ; dG%y_axis_units = oG%y_axis_units
   dG%x_ax_unit_short = oG%x_ax_unit_short ; dG%y_ax_unit_short = oG%y_ax_unit_short
+  dG%grid_unit_to_L = oG%grid_unit_to_L
   dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global
   dG%south_lat = oG%south_lat ; dG%west_lon  = oG%west_lon
   dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon

src/diagnostics/MOM_obsolete_params.F90

@@ -168,29 +168,15 @@ subroutine obsolete_logical(param_file, varname, warning_val, hint)
   character(len=*), optional, intent(in) :: hint   !< A hint to the user about what to do.
   ! Local variables
   logical :: test_logic, fatal_err
+  logical :: var_is_set  ! True if this value was read by read_param.
   character(len=128) :: hint_msg
 
-  test_logic = .false. ; call read_param(param_file, varname, test_logic)
+  test_logic = .false. ; call read_param(param_file, varname, test_logic, set=var_is_set)
   fatal_err = .true.
-  if (present(warning_val)) fatal_err = (warning_val .neqv. .true.)
+  if (var_is_set .and. present(warning_val)) fatal_err = (warning_val .neqv. test_logic)
   hint_msg = " " ; if (present(hint)) hint_msg = hint
 
-  if (test_logic) then
-    if (fatal_err) then
-      call MOM_ERROR(FATAL, "MOM_obsolete_params: "//trim(varname)//   &
-           " is an obsolete run-time flag, and should not be used. "// &
-           trim(hint_msg))
-    else
-      call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// &
-           " is an obsolete run-time flag. "//trim(hint_msg))
-    endif
-  endif
-
-  test_logic = .true. ; call read_param(param_file, varname, test_logic)
-  fatal_err = .true.
-  if (present(warning_val)) fatal_err = (warning_val .neqv. .false.)
-
-  if (.not.test_logic) then
+  if (var_is_set) then
     if (fatal_err) then
       call MOM_ERROR(FATAL, "MOM_obsolete_params: "//trim(varname)//   &
            " is an obsolete run-time flag, and should not be used. "// &
@@ -211,12 +197,13 @@ subroutine obsolete_char(param_file, varname, warning_val, hint)
   character(len=*), optional, intent(in) :: hint  !< A hint to the user about what to do.
   ! Local variables
   character(len=200) :: test_string, hint_msg
+  logical :: var_is_set  ! True if this value was read by read_param.
   logical :: only_warn
 
-  test_string = ''; call read_param(param_file, varname, test_string)
+  test_string = ''; call read_param(param_file, varname, test_string, set=var_is_set)
   hint_msg = " " ; if (present(hint)) hint_msg = hint
 
-  if (len_trim(test_string) > 0) then
+  if (var_is_set) then
     only_warn = .false.
     if (present(warning_val)) then ! Check if test_string and warning_val are the same.
       if (len_trim(warning_val) == len_trim(test_string)) then
@@ -246,15 +233,16 @@ subroutine obsolete_real(param_file, varname, warning_val, hint, only_warn)
 
   ! Local variables
   real :: test_val, warn_val
+  logical :: var_is_set  ! True if this value was read by read_param.
   logical :: issue_warning
   character(len=128) :: hint_msg
 
-  test_val = -9e35; call read_param(param_file, varname, test_val)
+  test_val = -9e35; call read_param(param_file, varname, test_val, set=var_is_set)
   warn_val = -9e35; if (present(warning_val)) warn_val = warning_val
   hint_msg = " " ; if (present(hint)) hint_msg = hint
   issue_warning = .false. ; if (present(only_warn)) issue_warning = only_warn
 
-  if (test_val /= -9e35) then
+  if (var_is_set) then
     if ((test_val == warn_val) .or. issue_warning) then
       call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// &
          " is an obsolete run-time flag. "//trim(hint_msg))
@@ -273,14 +261,15 @@ subroutine obsolete_int(param_file, varname, warning_val, hint)
   integer,     optional, intent(in) :: warning_val !< An allowed value that causes a warning instead of an error.
   character(len=*), optional, intent(in) :: hint   !< A hint to the user about what to do.
   ! Local variables
+  logical :: var_is_set  ! True if this value was read by read_param.
   integer :: test_val, warn_val
   character(len=128) :: hint_msg
 
-  test_val = -123456788; call read_param(param_file, varname, test_val)
+  test_val = -123456788; call read_param(param_file, varname, test_val, set=var_is_set)
   warn_val = -123456788; if (present(warning_val)) warn_val = warning_val
   hint_msg = " " ; if (present(hint)) hint_msg = hint
 
-  if (test_val /= -123456788) then
+  if (var_is_set) then
     if (test_val == warn_val) then
       call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// &
          " is an obsolete run-time flag. "//trim(hint_msg))

src/framework/MOM_dyn_horgrid.F90

@@ -181,6 +181,10 @@ module MOM_dyn_horgrid
 
   ! These parameters are run-time parameters that are used during some
   ! initialization routines (but not all)
+  real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related
+                        !! variables like len_lat and len_lon into rescaled horizontal distance
+                        !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or
+                        !! is 0 for a non-Cartesian grid.
   real :: south_lat     !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m]
   real :: west_lon      !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m]
   real :: len_lat       !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m]
@@ -400,6 +404,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns)
   G%len_lon = G_in%len_lon
 
   ! Rotation-invariant fields
+  G%grid_unit_to_L = G_in%grid_unit_to_L
   G%areaT_global = G_in%areaT_global
   G%IareaT_global = G_in%IareaT_global
   G%Rad_Earth = G_in%Rad_Earth