Merge branch 'dev/gfdl' into rescale_driver_time_vars

src/core/MOM_PressureForce_FV.F90

@@ -857,7 +857,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
   CS%GFS_scale = 1.0
   if (GV%g_prime(1) /= GV%g_Earth) CS%GFS_scale = GV%g_prime(1) / GV%g_Earth
 
-  call log_param(param_file, mdl, "GFS / G_EARTH", CS%GFS_scale)
+  call log_param(param_file, mdl, "GFS / G_EARTH", CS%GFS_scale, units="nondim")
 
 end subroutine PressureForce_FV_init
 

src/core/MOM_PressureForce_Montgomery.F90

@@ -874,7 +874,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_
   CS%GFS_scale = 1.0
   if (GV%g_prime(1) /= GV%g_Earth) CS%GFS_scale = GV%g_prime(1) / GV%g_Earth
 
-  call log_param(param_file, mdl, "GFS / G_EARTH", CS%GFS_scale)
+  call log_param(param_file, mdl, "GFS / G_EARTH", CS%GFS_scale, units="nondim")
 
 end subroutine PressureForce_Mont_init
 

src/core/MOM_barotropic.F90

@@ -4804,8 +4804,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
   endif
   if ((dtbt_tmp > 0.0) .and. (dtbt_input > 0.0)) calc_dtbt = .false.
 
-  call log_param(param_file, mdl, "DTBT as used", CS%dtbt*US%T_to_s)
-  call log_param(param_file, mdl, "estimated maximum DTBT", CS%dtbt_max*US%T_to_s)
+  call log_param(param_file, mdl, "DTBT as used", CS%dtbt*US%T_to_s, units="s")
+  call log_param(param_file, mdl, "estimated maximum DTBT", CS%dtbt_max*US%T_to_s, units="s")
 
   ! ubtav and vbtav, and perhaps ubt_IC and vbt_IC, are allocated and
   ! initialized in register_barotropic_restarts.

src/core/MOM_open_boundary.F90

@@ -5590,7 +5590,7 @@ end subroutine remap_OBC_fields
 !> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
 !!
 !! If the bottom most interface is below the topography then the bottom-most
-!! layers are contracted to GV%Angstrom_m.
+!! layers are contracted to GV%Angstrom_Z.
 !! If the bottom most interface is above the topography then the entire column
 !! is dilated (expanded) to fill the void.
 !!   @remark{There is a (hard-wired) "tolerance" parameter such that the

src/core/MOM_verticalGrid.F90

@@ -114,9 +114,9 @@ subroutine verticalGridInit( param_file, GV, US )
                  units="kg m-3", default=1035.0, scale=US%kg_m3_to_R)
   call get_param(param_file, mdl, "BOUSSINESQ", GV%Boussinesq, &
                  "If true, make the Boussinesq approximation.", default=.true.)
-  call get_param(param_file, mdl, "ANGSTROM", GV%Angstrom_m, &
+  call get_param(param_file, mdl, "ANGSTROM", GV%Angstrom_Z, &
                  "The minimum layer thickness, usually one-Angstrom.", &
-                 units="m", default=1.0e-10)
+                 units="m", default=1.0e-10, scale=US%m_to_Z)
   call get_param(param_file, mdl, "H_RESCALE_POWER", H_power, &
                  "An integer power of 2 that is used to rescale the model's "//&
                  "intenal units of thickness.  Valid values range from -300 to 300.", &
@@ -156,29 +156,29 @@ subroutine verticalGridInit( param_file, GV, US )
     GV%H_to_kg_m2 = US%R_to_kg_m3*GV%Rho0 * GV%H_to_m
     GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2
     GV%m_to_H = 1.0 / GV%H_to_m
-    GV%Angstrom_H = GV%m_to_H * GV%Angstrom_m
+    GV%Angstrom_H = GV%m_to_H * US%Z_to_m*GV%Angstrom_Z
     GV%H_to_MKS = GV%H_to_m
   else
     GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2
     GV%m_to_H = US%R_to_kg_m3*GV%Rho0 * GV%kg_m2_to_H
     GV%H_to_m = GV%H_to_kg_m2 / (US%R_to_kg_m3*GV%Rho0)
-    GV%Angstrom_H = GV%Angstrom_m*1000.0*GV%kg_m2_to_H
+    GV%Angstrom_H = US%Z_to_m*GV%Angstrom_Z * 1000.0*GV%kg_m2_to_H
     GV%H_to_MKS = GV%H_to_kg_m2
   endif
   GV%H_subroundoff = 1e-20 * max(GV%Angstrom_H,GV%m_to_H*1e-17)
   GV%H_to_Pa = US%L_T_to_m_s**2*US%m_to_Z * GV%g_Earth * GV%H_to_kg_m2
 
   GV%H_to_Z = GV%H_to_m * US%m_to_Z
   GV%Z_to_H = US%Z_to_m * GV%m_to_H
-  GV%Angstrom_Z = US%m_to_Z * GV%Angstrom_m
+  GV%Angstrom_m = US%Z_to_m * GV%Angstrom_Z
 
   GV%H_to_RZ = GV%H_to_kg_m2 * US%kg_m3_to_R * US%m_to_Z
   GV%RZ_to_H = GV%kg_m2_to_H * US%R_to_kg_m3 * US%Z_to_m
 
 ! Log derivative values.
-  call log_param(param_file, mdl, "M to THICKNESS", GV%m_to_H*H_rescale_factor)
-  call log_param(param_file, mdl, "M to THICKNESS rescaled by 2^-n", GV%m_to_H)
-  call log_param(param_file, mdl, "THICKNESS to M rescaled by 2^n", GV%H_to_m)
+  call log_param(param_file, mdl, "M to THICKNESS", GV%m_to_H*H_rescale_factor, units="H m-1")
+  call log_param(param_file, mdl, "M to THICKNESS rescaled by 2^-n", GV%m_to_H, units="2^n H m-1")
+  call log_param(param_file, mdl, "THICKNESS to M rescaled by 2^n", GV%H_to_m, units="2^-n m H-1")
 
   allocate( GV%sInterface(nk+1) )
   allocate( GV%sLayer(nk) )

src/framework/testing/MOM_file_parser_tests.F90

@@ -1277,7 +1277,7 @@ subroutine test_log_param_real
   call create_test_file(param_filename)
 
   call open_param_file(param_filename, param)
-  call log_param(param, module_name, sample_param_name, sample, desc=desc)
+  call log_param(param, module_name, sample_param_name, sample, desc=desc, units="")
   call close_param_file(param)
 end subroutine test_log_param_real
 
@@ -1290,7 +1290,7 @@ subroutine test_log_param_real_array
   call create_test_file(param_filename)
 
   call open_param_file(param_filename, param)
-  call log_param(param, module_name, sample_param_name, sample, desc=desc)
+  call log_param(param, module_name, sample_param_name, sample, desc=desc, units="")
   call close_param_file(param)
 end subroutine test_log_param_real_array
 

src/initialization/MOM_coord_initialization.F90

@@ -314,8 +314,8 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta
   real, dimension(GV%ke) :: T0   ! A profile of temperatures [C ~> degC]
   real, dimension(GV%ke) :: S0   ! A profile of salinities [S ~> ppt]
   real, dimension(GV%ke) :: Pref ! A array of reference pressures [R L2 T-2 ~> Pa]
-  real :: S_Ref   ! Default salinity range parameters [ppt].
-  real :: T_Ref   ! Default temperature range parameters [degC].
+  real :: S_Ref   ! Default salinity range parameters [S ~> ppt].
+  real :: T_Ref   ! Default temperature range parameters [C ~> degC].
   real :: S_Light, S_Dense ! Salinity range parameters [S ~> ppt].
   real :: T_Light, T_Dense ! Temperature range parameters [C ~> degC].
   real :: res_rat ! The ratio of density space resolution in the denser part
@@ -332,22 +332,26 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta
   call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
 
   call get_param(param_file, mdl, "T_REF", T_Ref, &
-                 "The default initial temperatures.", units="degC", default=10.0)
+                 "The default initial temperatures.", &
+                 units="degC", default=10.0, scale=US%degC_to_C)
   call get_param(param_file, mdl, "TS_RANGE_T_LIGHT", T_Light, &
                  "The initial temperature of the lightest layer when "//&
-                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref, scale=US%degC_to_C)
+                 "COORD_CONFIG is set to ts_range.", &
+                 units="degC", default=US%C_to_degC*T_Ref, scale=US%degC_to_C)
   call get_param(param_file, mdl, "TS_RANGE_T_DENSE", T_Dense, &
                  "The initial temperature of the densest layer when "//&
-                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref, scale=US%degC_to_C)
+                 "COORD_CONFIG is set to ts_range.", &
+                 units="degC", default=US%C_to_degC*T_Ref, scale=US%degC_to_C)
 
   call get_param(param_file, mdl, "S_REF", S_Ref, &
-                 "The default initial salinities.", units="PSU", default=35.0)
+                 "The default initial salinities.", &
+                 units="PSU", default=35.0, scale=US%ppt_to_S)
   call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", S_Light, &
-                 "The initial lightest salinities when COORD_CONFIG "//&
-                 "is set to ts_range.", default = S_Ref, units="PSU", scale=US%ppt_to_S)
+                 "The initial lightest salinities when COORD_CONFIG is set to ts_range.", &
+                 units="PSU", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S)
   call get_param(param_file, mdl, "TS_RANGE_S_DENSE", S_Dense, &
-                 "The initial densest salinities when COORD_CONFIG "//&
-                 "is set to ts_range.", default = S_Ref, units="PSU", scale=US%ppt_to_S)
+                 "The initial densest salinities when COORD_CONFIG is set to ts_range.", &
+                 units="PSU", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S)
 
   call get_param(param_file, mdl, "TS_RANGE_RESOLN_RATIO", res_rat, &
                  "The ratio of density space resolution in the densest "//&

src/initialization/MOM_state_initialization.F90

@@ -1782,8 +1782,8 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
                  "A reference temperature used in initialization.", &
                  units="degC", scale=US%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
   call get_param(param_file, mdl, "S_REF", S_Ref, &
-                 "A reference salinity used in initialization.", units="PSU", &
-                 default=35.0, scale=US%ppt_to_S, do_not_log=just_read)
+                 "A reference salinity used in initialization.", &
+                 units="PSU", default=35.0, scale=US%ppt_to_S, do_not_log=just_read)
   call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
                  "If true, accept the prescribed temperature and fit the "//&
                  "salinity; otherwise take salinity and fit temperature.", &
@@ -2480,8 +2480,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
   integer :: nkd      ! number of levels to use for regridding input arrays
   real    :: eps_Z    ! A negligibly thin layer thickness [Z ~> m].
   real    :: eps_rho  ! A negligibly small density difference [R ~> kg m-3].
-  real    :: PI_180   ! for conversion from degrees to radians
-  real    :: Hmix_default ! The default initial mixed layer depth [m].
+  real    :: PI_180   ! for conversion from degrees to radians [radian degree-1]
+  real    :: Hmix_default ! The default initial mixed layer depth [Z ~> m].
   real    :: Hmix_depth   ! The mixed layer depth in the initial condition [Z ~> m].
   real    :: missing_value_temp  ! The missing value in the input temperature field
   real    :: missing_value_salt  ! The missing value in the input salinity field
@@ -2680,10 +2680,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
                  default=.false., do_not_log=just_read.or.(GV%nkml==0))
     if (GV%nkml == 0) separate_mixed_layer = .false.
     call get_param(PF, mdl, "MINIMUM_DEPTH", Hmix_default, &
-                units="m", default=0.0, scale=1.0)
+                 units="m", default=0.0, scale=US%m_to_Z)
     call get_param(PF, mdl, "Z_INIT_HMIX_DEPTH", Hmix_depth, &
                  "The mixed layer depth in the initial conditions when Z_INIT_SEPARATE_MIXED_LAYER "//&
-                 "is set to true.", default=Hmix_default, units="m", scale=US%m_to_Z, &
+                 "is set to true.", units="m", default=US%Z_to_m*Hmix_default, scale=US%m_to_Z, &
                  do_not_log=(just_read .or. .not.separate_mixed_layer))
     call get_param(PF, mdl, "LAYER_Z_INIT_IC_EXTRAP_BUG", density_extrap_bug, &
                  "If true use an expression with a vertical indexing bug for extrapolating the "//&

src/parameterizations/lateral/MOM_MEKE.F90

@@ -1196,8 +1196,8 @@ logical function MEKE_init(Time, G, US, param_file, diag, dbcomms_CS, CS, MEKE,
                    "each time step.", default=.false.)
     if (CS%MEKE_equilibrium_restoring) then
       call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, &
-                     "The timescale used to nudge MEKE toward its equilibrium value.", units="s", &
-                     default=1e6, scale=US%s_to_T)
+                     "The timescale used to nudge MEKE toward its equilibrium value.", &
+                     units="s", default=1e6, scale=US%s_to_T)
       CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale
     endif
 
@@ -1210,8 +1210,8 @@ logical function MEKE_init(Time, G, US, param_file, diag, dbcomms_CS, CS, MEKE,
                    "by GME.  If MEKE_GMECOEFF is negative, this conversion "//&
                    "is not used or calculated.", units="nondim", default=-1.0)
     call get_param(param_file, mdl, "MEKE_BGSRC", CS%MEKE_BGsrc, &
-                   "A background energy source for MEKE.", units="W kg-1", &
-                   default=0.0, scale=US%m_to_L**2*US%T_to_s**3)
+                   "A background energy source for MEKE.", &
+                   units="W kg-1", default=0.0, scale=US%m_to_L**2*US%T_to_s**3)
     call get_param(param_file, mdl, "MEKE_KH", CS%MEKE_Kh, &
                    "A background lateral diffusivity of MEKE. "//&
                    "Use a negative value to not apply lateral diffusion to MEKE.", &
@@ -1248,11 +1248,9 @@ logical function MEKE_init(Time, G, US, param_file, diag, dbcomms_CS, CS, MEKE,
                  "If true, use the vertvisc_type to calculate the bottom "//&
                  "drag acting on MEKE.", default=.true.)
   call get_param(param_file, mdl, "MEKE_KHTH_FAC", MEKE%KhTh_fac, &
-                 "A factor that maps MEKE%Kh to KhTh.", units="nondim", &
-                 default=0.0)
+                 "A factor that maps MEKE%Kh to KhTh.", units="nondim", default=0.0)
   call get_param(param_file, mdl, "MEKE_KHTR_FAC", MEKE%KhTr_fac, &
-                 "A factor that maps MEKE%Kh to KhTr.", units="nondim", &
-                 default=0.0)
+                 "A factor that maps MEKE%Kh to KhTr.", units="nondim", default=0.0)
   call get_param(param_file, mdl, "MEKE_KHMEKE_FAC", CS%KhMEKE_Fac, &
                  "A factor that maps MEKE%Kh to Kh for MEKE itself.", &
                  units="nondim", default=0.0)
@@ -1336,13 +1334,11 @@ logical function MEKE_init(Time, G, US, param_file, diag, dbcomms_CS, CS, MEKE,
 
   ! Nonlocal module parameters
   call get_param(param_file, mdl, "CDRAG", cdrag, &
-                 "CDRAG is the drag coefficient relating the magnitude of "//&
-                 "the velocity field to the bottom stress.", units="nondim", &
-                 default=0.003)
+                 "CDRAG is the drag coefficient relating the magnitude of the velocity "//&
+                 "field to the bottom stress.", units="nondim", default=0.003)
   call get_param(param_file, mdl, "MEKE_CDRAG", CS%cdrag, &
                  "Drag coefficient relating the magnitude of the velocity "//&
-                 "field to the bottom stress in MEKE.", units="nondim", &
-                 default=cdrag)
+                 "field to the bottom stress in MEKE.", units="nondim", default=cdrag)
   call get_param(param_file, mdl, "LAPLACIAN", laplacian, default=.false., do_not_log=.true.)
   call get_param(param_file, mdl, "BIHARMONIC", biharmonic, default=.false., do_not_log=.true.)
 

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -1194,7 +1194,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
                  default=1.0e-17, units="s-1", scale=US%T_to_s)
   call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_FGNV_streamfn, &
                  default=.false., do_not_log=.true.)
-  CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn
+  CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct
   CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE
   ! Indicate whether to calculate the Eady growth rate
   CS%calculate_Eady_growth_rate = use_MEKE .or. (KhTr_Slope_Cff>0.) .or. (KhTh_Slope_Cff>0.)

src/parameterizations/lateral/MOM_tidal_forcing.F90

@@ -70,16 +70,16 @@ module MOM_tidal_forcing
     sin_struct(:,:,:), &    !< The sine and cosine based structures that can
     cos_struct(:,:,:), &    !< be associated with the astronomical forcing [nondim].
     cosphasesal(:,:,:), &   !< The cosine and sine of the phase of the
-    sinphasesal(:,:,:), &   !< self-attraction and loading amphidromes.
+    sinphasesal(:,:,:), &   !< self-attraction and loading amphidromes [nondim].
     ampsal(:,:,:), &        !< The amplitude of the SAL [Z ~> m].
     cosphase_prev(:,:,:), & !< The cosine and sine of the phase of the
-    sinphase_prev(:,:,:), & !< amphidromes in the previous tidal solutions.
+    sinphase_prev(:,:,:), & !< amphidromes in the previous tidal solutions [nondim].
     amp_prev(:,:,:)         !< The amplitude of the previous tidal solution [Z ~> m].
   type(sht_CS) :: sht       !< Spherical harmonic transforms (SHT) for SAL
-  integer :: sal_sht_Nd     !< Maximum degree for SHT [nodim]
-  real, allocatable :: Love_Scaling(:)      !< Love number for each SHT mode [nodim]
-  real, allocatable :: Snm_Re(:), & !< Real and imaginary SHT coefficient for SHT SAL
-                       Snm_Im(:)    !< [Z ~> m]
+  integer :: sal_sht_Nd     !< Maximum degree for SHT [nondim]
+  real, allocatable :: Love_Scaling(:)  !< Love number for each SHT mode [nondim]
+  real, allocatable :: Snm_Re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m]
+                       Snm_Im(:)    !< Imaginary SHT coefficient for SHT SAL [Z ~> m]
 end type tidal_forcing_CS
 
 integer :: id_clock_tides !< CPU clock for tides
@@ -99,9 +99,12 @@ module MOM_tidal_forcing
 subroutine astro_longitudes_init(time_ref, longitudes)
   type(time_type), intent(in) :: time_ref            !> Time to calculate longitudes for.
   type(astro_longitudes), intent(out) :: longitudes  !> Lunar and solar longitudes at time_ref.
+
+  ! Local variables
   real :: D                                          !> Time since the reference date [days]
   real :: T                                          !> Time in Julian centuries [centuries]
   real, parameter :: PI = 4.0 * atan(1.0)            !> 3.14159... [nondim]
+
   ! Find date at time_ref in days since 1900-01-01
   D = time_type_to_real(time_ref - set_date(1900, 1, 1)) / (24.0 * 3600.0)
   ! Time since 1900-01-01 in Julian centuries
@@ -542,7 +545,7 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS)
     call get_param(param_file, mdl, "TIDAL_SAL_SHT_DEGREE", CS%sal_sht_Nd, &
                    "The maximum degree of the spherical harmonics transformation used for "// &
                    "calculating the self-attraction and loading term for tides.", &
-                   default=0, do_not_log=.not. CS%tidal_sal_sht)
+                   default=0, do_not_log=.not.CS%tidal_sal_sht)
     call get_param(param_file, mdl, "RHO_0", rhoW, &
                    "The mean ocean density used with BOUSSINESQ true to "//&
                    "calculate accelerations and the mass for conservation "//&
@@ -551,8 +554,9 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS)
                    units="kg m-3", default=1035.0, scale=US%kg_m3_to_R, do_not_log=.True.)
     call get_param(param_file, mdl, "RHO_E", rhoE, &
                    "The mean solid earth density.  This is used for calculating the "// &
-                   "self-attraction and loading term.", units="kg m-3", &
-                   default=5517.0, scale=US%kg_m3_to_R, do_not_log=.not. CS%tidal_sal_sht)
+                   "self-attraction and loading term.", &
+                   units="kg m-3", default=5517.0, scale=US%kg_m3_to_R, &
+                   do_not_log=.not.CS%tidal_sal_sht)
     lmax = calc_lmax(CS%sal_sht_Nd)
     allocate(CS%Snm_Re(lmax)); CS%Snm_Re(:) = 0.0
     allocate(CS%Snm_Im(lmax)); CS%Snm_Im(:) = 0.0