Merge branch 'dev/gfdl' into rescale_face_lengths_list

config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90

@@ -28,6 +28,7 @@ module g_tracer_utils
     character(len=fm_string_len) :: obc_src_file_name !< Boundary condition tracer source filename
     character(len=fm_string_len) :: obc_src_field_name !< Boundary condition tracer source fieldname
     integer :: src_var_record !< Unknown
+    logical :: runoff_added_to_stf = .false. !< Has flux in from runoff been added to stf?
     logical :: requires_src_info = .false. !< Unknown
     real    :: src_var_unit_conversion = 1.0 !< This factor depends on the tracer. Ask Jasmin
     real    :: src_var_valid_min = 0.0 !< Unknown

src/ALE/MOM_regridding.F90

@@ -179,6 +179,9 @@ module MOM_regridding
 !> Default minimum thickness for some coordinate generation modes
 real, parameter, public :: regriddingDefaultMinThickness = 1.e-3
 
+!> Maximum length of parameters
+integer, parameter :: MAX_PARAM_LENGTH = 120
+
 #undef __DO_SAFETY_CHECKS__
 
 contains
@@ -199,7 +202,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
   ! Local variables
   integer :: ke ! Number of levels
   character(len=80)  :: string, string2, varName ! Temporary strings
-  character(len=40)  :: coord_units, param_name, coord_res_param ! Temporary strings
+  character(len=40)  :: coord_units, coord_res_param ! Temporary strings
+  character(len=MAX_PARAM_LENGTH) :: param_name
   character(len=200) :: inputdir, fileName
   character(len=320) :: message ! Temporary strings
   character(len=12) :: expected_units, alt_units ! Temporary strings
@@ -256,7 +260,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
       param_name = "INTERPOLATION_SCHEME"
       string2 = regriddingDefaultInterpScheme
     else
-      param_name = trim(param_prefix)//"_INTERP_SCHEME_"//trim(param_suffix)
+      param_name = create_coord_param(param_prefix, "INTERP_SCHEME", param_suffix)
       string2 = 'PPM_H4' ! Default for diagnostics
     endif
     call get_param(param_file, mdl, "INTERPOLATION_SCHEME", string, &
@@ -309,8 +313,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
     coord_res_param = "ALE_RESOLUTION"
     string2 = 'UNIFORM'
   else
-    param_name = trim(param_prefix)//"_DEF_"//trim(param_suffix)
-    coord_res_param = trim(param_prefix)//"_RES_"//trim(param_suffix)
+    param_name = create_coord_param(param_prefix, "DEF", param_suffix)
+    coord_res_param = create_coord_param(param_prefix, "RES", param_suffix)
     string2 = 'UNIFORM'
     if (maximum_depth>3000.) string2='WOA09' ! For convenience
   endif
@@ -545,13 +549,22 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
   ! initialise coordinate-specific control structure
   call initCoord(CS, GV, US, coord_mode, param_file)
 
-  if (main_parameters .and. coord_is_state_dependent) then
-    call get_param(param_file, mdl, "P_REF", P_Ref, &
-                 "The pressure that is used for calculating the coordinate "//&
-                 "density.  (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
-                 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
-                 units="Pa", default=2.0e7, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
-    call get_param(param_file, mdl, "REGRID_COMPRESSIBILITY_FRACTION", tmpReal, &
+  if (coord_is_state_dependent) then
+    if (main_parameters) then
+      call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), P_Ref, &
+                   "The pressure that is used for calculating the coordinate "//&
+                   "density.  (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
+                   "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
+                   units="Pa", default=2.0e7, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
+    else
+      call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), P_Ref, &
+                   "The pressure that is used for calculating the diagnostic coordinate "//&
+                   "density.  (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
+                   "This is only used for the RHO coordinate.", &
+                   units="Pa", default=2.0e7, scale=US%kg_m3_to_R*US%m_s_to_L_T**2)
+    endif
+    call get_param(param_file, mdl, create_coord_param(param_prefix, "REGRID_COMPRESSIBILITY_FRACTION", param_suffix), &
+                 tmpReal, &
                  "When interpolating potential density profiles we can add "//&
                  "some artificial compressibility solely to make homogeneous "//&
                  "regions appear stratified.", units="nondim", default=0.)
@@ -2432,6 +2445,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri
     if (present(min_thickness)) call set_sigma_params(CS%sigma_CS, min_thickness=min_thickness)
   case (REGRIDDING_RHO)
     if (present(min_thickness)) call set_rho_params(CS%rho_CS, min_thickness=min_thickness)
+    if (present(ref_pressure)) call set_rho_params(CS%rho_CS, ref_pressure=ref_pressure)
     if (present(integrate_downward_for_e)) &
       call set_rho_params(CS%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
     if (associated(CS%rho_CS) .and. (present(interp_scheme) .or. present(boundary_extrapolation))) &
@@ -2564,6 +2578,29 @@ subroutine dz_function1( string, dz )
 
 end subroutine dz_function1
 
+!> Construct the name of a parameter for a specific coordinate based on param_prefix and param_suffix. For the main,
+!! prognostic coordinate this will simply return the parameter name (e.g. P_REF)
+function create_coord_param(param_prefix, param_name, param_suffix) result(coord_param)
+  character(len=*) :: param_name   !< The base name of the parameter (e.g. the one used for the main coordinate)
+  character(len=*) :: param_prefix !< String to prefix to parameter names.
+  character(len=*) :: param_suffix !< String to append to parameter names.
+  character(len=MAX_PARAM_LENGTH) :: coord_param  !< Parameter name prepended by param_prefix
+                                                  !! and appended with param_suffix
+  integer :: out_length
+
+  if (len_trim(param_prefix) + len_trim(param_suffix) == 0) then
+    coord_param = param_name
+  else
+    ! Note the +2 is because of two underscores
+    out_length = len_trim(param_name)+len_trim(param_prefix)+len_trim(param_suffix)+2
+    if (out_length > MAX_PARAM_LENGTH) then
+      call MOM_error(FATAL,"Coordinate parameter is too long; increase MAX_PARAM_LENGTH")
+    endif
+    coord_param = TRIM(param_prefix)//"_"//TRIM(param_name)//"_"//TRIM(param_suffix)
+  endif
+
+end function create_coord_param
+
 !> Parses a string and generates a rho_target(:) profile with refined resolution downward
 !! and returns the number of levels
 integer function rho_function1( string, rho_target )

src/ALE/coord_rho.F90

@@ -67,12 +67,14 @@ subroutine end_coord_rho(CS)
 end subroutine end_coord_rho
 
 !> This subroutine can be used to set the parameters for the coord_rho module
-subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS)
+subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS, ref_pressure)
   type(rho_CS),      pointer    :: CS !< Coordinate control structure
   real,    optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
   logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface
                                       !! positions from the top downward.  If false, integrate
                                       !! from the bottom upward, as does the rest of the model.
+  real,    optional, intent(in) :: ref_pressure     !< The reference pressure for density-dependent
+                                                    !! coordinates [R L2 T-2 ~> Pa]
 
   type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation
 
@@ -81,6 +83,7 @@ subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS
   if (present(min_thickness)) CS%min_thickness = min_thickness
   if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e
   if (present(interp_CS)) CS%interp_CS = interp_CS
+  if (present(ref_pressure)) CS%ref_pressure = ref_pressure
 end subroutine set_rho_params
 
 !> Build a rho coordinate column

src/core/MOM.F90

@@ -114,7 +114,8 @@ module MOM
 use MOM_shared_initialization, only : write_ocean_geometry_file
 use MOM_sponge,                only : init_sponge_diags, sponge_CS
 use MOM_state_initialization,  only : MOM_initialize_state
-use MOM_stoch_eos,             only : MOM_stoch_eos_init,MOM_stoch_eos_run,MOM_stoch_eos_CS,mom_calc_varT
+use MOM_stoch_eos,             only : MOM_stoch_eos_init, MOM_stoch_eos_run, MOM_stoch_eos_CS
+use MOM_stoch_eos,             only : stoch_EOS_register_restarts, post_stoch_EOS_diags, mom_calc_varT
 use MOM_sum_output,            only : write_energy, accumulate_net_input
 use MOM_sum_output,            only : MOM_sum_output_init, MOM_sum_output_end
 use MOM_sum_output,            only : sum_output_CS
@@ -288,6 +289,7 @@ module MOM
   logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics.
   logical :: mixedlayer_restrat      !< If true, use submesoscale mixed layer restratifying scheme.
   logical :: useMEKE                 !< If true, call the MEKE parameterization.
+  logical :: use_stochastic_EOS      !< If true, use the stochastic EOS parameterizations.
   logical :: useWaves                !< If true, update Stokes drift
   logical :: use_p_surf_in_EOS       !< If true, always include the surface pressure contributions
                                      !! in equation of state calculations.
@@ -298,7 +300,10 @@ module MOM
                                      !! calculated, and if it is 0, dtbt is calculated every step.
   type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period.
   type(time_type) :: dtbt_reset_time     !< The next time DTBT should be calculated.
-  real            :: dt_obc_seg_period   !< The time interval between OBC segment updates for OBGC tracers
+  real            :: dt_obc_seg_period   !< The time interval between OBC segment updates for OBGC
+                                         !! tracers [T ~> s], or a negative value if the segment
+                                         !! data are time-invarant, or zero to update the OBGC
+                                         !! segment data with every call to update_OBC_segment_data.
   type(time_type) :: dt_obc_seg_interval !< A time_time representation of dt_obc_seg_period.
   type(time_type) :: dt_obc_seg_time     !< The next time OBC segment update is applied to OBGC tracers.
 
@@ -1079,12 +1084,12 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
 
   call cpu_clock_begin(id_clock_dynamics)
   call cpu_clock_begin(id_clock_stoch)
-  if (CS%stoch_eos_CS%use_stoch_eos) call MOM_stoch_eos_run(G,u,v,dt,Time_local,CS%stoch_eos_CS,CS%diag)
+  if (CS%use_stochastic_EOS) call MOM_stoch_eos_run(G, u, v, dt, Time_local, CS%stoch_eos_CS)
   call cpu_clock_end(id_clock_stoch)
   call cpu_clock_begin(id_clock_varT)
-  if (CS%stoch_eos_CS%stanley_coeff >= 0.0) then
-    call MOM_calc_varT(G,GV,h,CS%tv,CS%stoch_eos_CS,dt)
-    call pass_var(CS%tv%varT, G%Domain,clock=id_clock_pass,halo=1)
+  if (CS%use_stochastic_EOS) then
+    call MOM_calc_varT(G, GV, h, CS%tv, CS%stoch_eos_CS, dt)
+    if (associated(CS%tv%varT)) call pass_var(CS%tv%varT, G%Domain, clock=id_clock_pass, halo=1)
   endif
   call cpu_clock_end(id_clock_varT)
 
@@ -1297,9 +1302,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   if (IDs%id_u > 0) call post_data(IDs%id_u, u, CS%diag)
   if (IDs%id_v > 0) call post_data(IDs%id_v, v, CS%diag)
   if (IDs%id_h > 0) call post_data(IDs%id_h, h, CS%diag)
-  if (CS%stoch_eos_CS%id_stoch_eos > 0) call post_data(CS%stoch_eos_CS%id_stoch_eos, CS%stoch_eos_CS%pattern, CS%diag)
-  if (CS%stoch_eos_CS%id_stoch_phi > 0) call post_data(CS%stoch_eos_CS%id_stoch_phi, CS%stoch_eos_CS%phi, CS%diag)
-  if (CS%stoch_eos_CS%id_tvar_sgs > 0) call post_data(CS%stoch_eos_CS%id_tvar_sgs, CS%tv%varT, CS%diag)
+  if (CS%use_stochastic_EOS) call post_stoch_EOS_diags(CS%stoch_eos_CS, CS%tv, CS%diag)
   call disable_averaging(CS%diag)
   call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
 
@@ -2186,12 +2189,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                  units="s", default=default_val, do_not_read=(dtbt > 0.0))
   endif
 
-  CS%dt_obc_seg_period = -1.0
   call get_param(param_file, "MOM", "DT_OBC_SEG_UPDATE_OBGC", CS%dt_obc_seg_period, &
                "The time between OBC segment data updates for OBGC tracers. "//&
                "This must be an integer multiple of DT and DT_THERM. "//&
                "The default is set to DT.", &
-               units="s", default=US%T_to_s*CS%dt, do_not_log=.not.associated(CS%OBC))
+               units="s", default=US%T_to_s*CS%dt, scale=US%s_to_T, do_not_log=.not.associated(CS%OBC))
 
   ! This is here in case these values are used inappropriately.
   use_frazil = .false. ; bound_salinity = .false.
@@ -2219,11 +2221,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                  "A tiny magnitude of temperatures below which they are set to 0.", &
                  units="degC", default=0.0, scale=US%degC_to_C)
     call get_param(param_file, "MOM", "C_P", CS%tv%C_p, &
-                 "The heat capacity of sea water, approximated as a "//&
-                 "constant. This is only used if ENABLE_THERMODYNAMICS is "//&
-                 "true. The default value is from the TEOS-10 definition "//&
-                 "of conservative temperature.", units="J kg-1 K-1", &
-                 default=3991.86795711963, scale=US%J_kg_to_Q*US%C_to_degC)
+                 "The heat capacity of sea water, approximated as a constant. "//&
+                 "This is only used if ENABLE_THERMODYNAMICS is true. The default "//&
+                 "value is from the TEOS-10 definition of conservative temperature.", &
+                 units="J kg-1 K-1", default=3991.86795711963, scale=US%J_kg_to_Q*US%C_to_degC)
     call get_param(param_file, "MOM", "USE_PSURF_IN_EOS", CS%use_p_surf_in_EOS, &
                  "If true, always include the surface pressure contributions "//&
                  "in equation of state calculations.", default=.true.)
@@ -2239,9 +2240,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                  "The number of sublayers within the mixed layer if "//&
                  "BULKMIXEDLAYER is true.", units="nondim", default=2)
     call get_param(param_file, "MOM", "NKBL", nkbl, &
-                 "The number of layers that are used as variable density "//&
-                 "buffer layers if BULKMIXEDLAYER is true.", units="nondim", &
-                 default=2)
+                 "The number of layers that are used as variable density buffer "//&
+                 "layers if BULKMIXEDLAYER is true.", units="nondim", default=2)
   endif
 
   call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, &
@@ -2642,7 +2642,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
 
   call MEKE_alloc_register_restart(HI, US, param_file, CS%MEKE, restart_CSp)
   call set_visc_register_restarts(HI, GV, US, param_file, CS%visc, restart_CSp)
-  call mixedlayer_restrat_register_restarts(HI, GV, param_file, &
+  call mixedlayer_restrat_register_restarts(HI, GV, US, param_file, &
            CS%mixedlayer_restrat_CSp, restart_CSp)
 
   if (CS%rotate_index .and. associated(OBC_in) .and. use_temperature) then
@@ -2678,7 +2678,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   endif
 
   if (present(waves_CSp)) then
-    call waves_register_restarts(waves_CSp, HI, GV, param_file, restart_CSp)
+    call waves_register_restarts(waves_CSp, HI, GV, US, param_file, restart_CSp)
+  endif
+
+  if (use_temperature) then
+    call stoch_EOS_register_restarts(HI, param_file, CS%stoch_eos_CS, restart_CSp)
   endif
 
   call callTree_waypoint("restart registration complete (initialize_MOM)")
@@ -2966,7 +2970,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     call interface_filter_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%interface_filter_CSp)
 
   new_sim = is_new_run(restart_CSp)
-  call MOM_stoch_eos_init(G,Time,param_file,CS%stoch_eos_CS,restart_CSp,diag)
+  if (use_temperature) then
+    CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, US, param_file, diag, CS%stoch_eos_CS, restart_CSp)
+  else
+    CS%use_stochastic_EOS = .false.
+  endif
 
   if (CS%use_porbar) &
     call porous_barriers_init(Time, US, param_file, diag, CS%por_bar_CS)
@@ -3209,7 +3217,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp)
   type(unit_scale_type),   pointer :: US => NULL() ! Pointer to a structure containing
                                                    ! various unit conversion factors
   type(MOM_restart_CS),    pointer :: restart_CSp_tmp => NULL()
-  real, allocatable :: z_interface(:,:,:) ! Interface heights [m]
+  real, allocatable :: z_interface(:,:,:) ! Interface heights [Z ~> m]
 
   call cpu_clock_begin(id_clock_init)
   call callTree_enter("finish_MOM_initialization()")
@@ -3232,9 +3240,9 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp)
     restart_CSp_tmp = restart_CSp
     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, eta_to_m=1.0, dZref=G%Z_ref)
+    call find_eta(CS%h, CS%tv, G, GV, US, z_interface, dZref=G%Z_ref)
     call register_restart_field(z_interface, "eta", .true., restart_CSp_tmp, &
-                                "Interface heights", "meter", z_grid='i')
+                                "Interface heights", "meter", z_grid='i', conversion=US%Z_to_m)
     ! NOTE: write_ic=.true. routes routine to fms2 IO write_initial_conditions interface
     call save_restart(dirs%output_directory, Time, CS%G_in, &
                       restart_CSp_tmp, filename=CS%IC_file, GV=GV, write_ic=.true.)