From dcdc509c1c91769e68e944830cc84aaa7ef892ac Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 14 Nov 2018 17:46:07 -0500 Subject: [PATCH 01/11] +Add US arg to set_up_ALE_sponge_vel_field_varying Use US%m_to_Z in place of G%Zd_to_m to convert units from m to Z in set_up_ALE_sponge_vel_field_varying, whic required adding a unit_scale_type argument. All answers are bitwise identical. --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 1ad8e44960..b08e77e213 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -741,13 +741,14 @@ end subroutine set_up_ALE_sponge_vel_field_fixed !> This subroutine stores the reference profile at uand v points for the variable !! whose address is given by u_ptr and v_ptr. subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename_v, fieldname_v, & - Time, G, CS, u_ptr, v_ptr) + Time, G, US, CS, u_ptr, v_ptr) character(len=*), intent(in) :: filename_u !< File name for u field character(len=*), intent(in) :: fieldname_u !< Name of u variable in file character(len=*), intent(in) :: filename_v !< File name for v field character(len=*), intent(in) :: fieldname_v !< Name of v variable in file type(time_type), intent(in) :: Time !< Model time type(ocean_grid_type), intent(inout) :: G !< Ocean grid (in) + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(in) :: u_ptr !< u pointer to the field to be damped (in). real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(in) :: v_ptr !< v pointer to the field to be damped (in). @@ -798,7 +799,7 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename ! modulo attribute of the zonal axis (mjh). call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id,Time, 1.0,G,u_val,mask_u,z_in,z_edges_in,& - missing_value,.true.,.false.,.false., m_to_Z=1.0/G%Zd_to_m) + missing_value,.true.,.false.,.false., m_to_Z=US%m_to_Z) !!! TODO: add a velocity interface! (mjh) @@ -808,7 +809,7 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename ! modulo attribute of the zonal axis (mjh). call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id,Time, 1.0,G,v_val,mask_v,z_in,z_edges_in, & - missing_value,.true.,.false.,.false., m_to_Z=1.0/G%Zd_to_m) + missing_value,.true.,.false.,.false., m_to_Z=US%m_to_Z) ! stores the reference profile allocate(CS%Ref_val_u%p(fld_sz(3),CS%num_col_u)) From 375bab5a0ea0ae662e100256ad5c26e95ddd1dbf Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 14 Nov 2018 17:46:57 -0500 Subject: [PATCH 02/11] +Rescaled variables in MOM_internal_tides Recast internal depth and height variables in MOM_internal_tides to use units of Z instead of m. This required adding a unit_scale_type arguments to internal_tides_init and itidal_lowmode_loss, and a modified call from diabatic_driver_init. All answers are bitwise identical. --- .../lateral/MOM_internal_tides.F90 | 22 +++++++++++-------- .../vertical/MOM_diabatic_driver.F90 | 2 +- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 8978a9b1f9..527eca30bc 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -74,7 +74,7 @@ module MOM_internal_tides !< energy lost due to wave breaking [W m-2] real, allocatable, dimension(:,:) :: TKE_itidal_loss_fixed !< fixed part of the energy lost due to small-scale drag - !! [kg m-2] here; will be multiplied by N and En to get into [W m-2] + !! [kg Z-2] here; will be multiplied by N and En to get into [W m-2] real, allocatable, dimension(:,:,:,:,:) :: TKE_itidal_loss !< energy lost due to small-scale wave drag [W m-2] real, allocatable, dimension(:,:) :: tot_leak_loss !< Energy loss rates due to misc bakground processes, @@ -397,7 +397,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! Finally, apply loss if (CS%apply_wave_drag) then ! Calculate loss rate and apply loss over the time step - call itidal_lowmode_loss(G, CS, Nb, Ub, CS%En, CS%TKE_itidal_loss_fixed, & + call itidal_lowmode_loss(G, US, CS, Nb, Ub, CS%En, CS%TKE_itidal_loss_fixed, & CS%TKE_itidal_loss, dt, full_halos=.false.) endif ! Check for En<0 - for debugging, delete later @@ -622,8 +622,9 @@ end subroutine sum_En !> Calculates the energy lost from the propagating internal tide due to !! scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001). -subroutine itidal_lowmode_loss(G, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos) +subroutine itidal_lowmode_loss(G, US, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(int_tide_CS), pointer :: CS !< The control structure returned by a !! previous call to int_tide_init. real, dimension(G%isd:G%ied,G%jsd:G%jed), & @@ -633,7 +634,7 @@ subroutine itidal_lowmode_loss(G, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, !! mode velocity, in m s-1. real, dimension(G%isd:G%ied,G%jsd:G%jed), & intent(in) :: TKE_loss_fixed !< Fixed part of energy loss, - !! in kg m-2 (rho*kappa*h^2). + !! in kg Z-2 (rho*kappa*h^2). real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), & intent(inout) :: En !< Energy density of the internal waves, in J m-2. real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), & @@ -670,7 +671,7 @@ subroutine itidal_lowmode_loss(G, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, enddo ! Calculate TKE loss rate; units of [W m-2] here. - TKE_loss_tot = q_itides * TKE_loss_fixed(i,j) * Nb(i,j) * Ub(i,j,fr,m)**2 + TKE_loss_tot = q_itides * US%Z_to_m**2 * TKE_loss_fixed(i,j) * Nb(i,j) * Ub(i,j,fr,m)**2 ! Update energy remaining (this is a pseudo implicit calc) ! (E(t+1)-E(t))/dt = -TKE_loss(E(t+1)/E(t)), which goes to zero as E(t+1) goes to zero @@ -2104,10 +2105,11 @@ end subroutine PPM_limit_pos ! end subroutine register_int_tide_restarts !> This subroutine initializes the internal tides module. -subroutine internal_tides_init(Time, G, GV, param_file, diag, CS) +subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) type(time_type), target, intent(in) :: Time !< The current model time. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time !! parameters. type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate @@ -2280,7 +2282,7 @@ subroutine internal_tides_init(Time, G, GV, param_file, diag, CS) "dissipated locally with INT_TIDE_DISSIPATION. \n"//& "THIS NAME COULD BE BETTER.", & units="nondim", default=0.3333) - call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, & + call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, & "A topographic wavenumber used with INT_TIDE_DISSIPATION. \n"//& "The default is 2pi/10 km, as in St.Laurent et al. 2002.", & units="m-1", default=8.e-4*atan(1.0)) @@ -2312,16 +2314,18 @@ subroutine internal_tides_init(Time, G, GV, param_file, diag, CS) fail_if_missing=.true.) filename = trim(CS%inputdir) // trim(h2_file) call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename) - call MOM_read_data(filename, 'h2', h2, G%domain, timelevel=1) + call MOM_read_data(filename, 'h2', h2, G%domain, timelevel=1, scale=US%m_to_Z) do j=G%jsc,G%jec ; do i=G%isc,G%iec ! Restrict rms topo to 10 percent of column depth. - h2(i,j) = min(0.01*(G%Zd_to_m*G%bathyT(i,j))**2, h2(i,j)) + h2(i,j) = min(0.01*(G%bathyT(i,j))**2, h2(i,j)) ! Compute the fixed part; units are [kg m-2] here ! will be multiplied by N and En to get into [W m-2] CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*GV%Rho0*& kappa_itides * h2(i,j) enddo ; enddo + deallocate(h2) + ! Read in prescribed coast/ridge/shelf angles from file call get_param(param_file, mdl, "REFL_ANGLE_FILE", refl_angle_file, & "The path to the file containing the local angle of \n"//& diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 981fcd0b6a..4506177f41 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3295,7 +3295,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di if (CS%use_int_tides) then call int_tide_input_init(Time, G, GV, US, param_file, diag, CS%int_tide_input_CSp, & CS%int_tide_input) - call internal_tides_init(Time, G, GV, param_file, diag, CS%int_tide_CSp) + call internal_tides_init(Time, G, GV, US, param_file, diag, CS%int_tide_CSp) endif ! initialize module for setting diffusivities From 967e470ebe6f9b332f7133ac5a95158d3c0afeb0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 14 Nov 2018 17:48:04 -0500 Subject: [PATCH 03/11] +Rescaled variables in MOM_tidal_mixing Recast internal depth and height variables in MOM_tidal_mixing to use units of Z instead of m. This required adding a unit_scale_type arguments to the internal subroutines read_tidal_energy and read_tidal_constituents. Also eliminated the use of where statements and array syntax, which are strongly discouraged because they can operate on uninitialized points in the halo regions. All answers are bitwise identical. --- .../vertical/MOM_tidal_mixing.F90 | 92 ++++++++++--------- 1 file changed, 49 insertions(+), 43 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 254843ebe2..05c377ab9c 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -136,8 +136,8 @@ module MOM_tidal_mixing type(CVMix_tidal_params_type) :: CVMix_tidal_params !< A CVMix-specific type with parameters for tidal mixing type(CVMix_global_params_type) :: CVMix_glb_params !< CVMix-specific for Prandtl number only real :: tidal_max_coef !< CVMix-specific maximum allowable tidal diffusivity. [m^2/s] - real :: tidal_diss_lim_tc !< CVMix-specific dissipation limit for - !! tidal-energy-constituent data + real :: tidal_diss_lim_tc !< CVMix-specific dissipation limit depth for + !! tidal-energy-constituent data, in Z. type(remapping_CS) :: remap_CS !< The control structure for remapping ! Data containers @@ -510,7 +510,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ call get_param(param_file, mdl, "TIDAL_DISS_LIM_TC", CS%tidal_diss_lim_tc, & "Min allowable depth for dissipation for tidal-energy-constituent data. \n"//& "No dissipation contribution is applied above TIDAL_DISS_LIM_TC.", & - units="m", default=0.0) + units="m", default=0.0, scale=US%m_to_Z) call get_param(param_file, mdl, "TIDAL_ENERGY_FILE",tidal_energy_file, & "The path to the file containing tidal energy \n"//& "dissipation. Used with CVMix tidal mixing schemes.", & @@ -549,7 +549,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, diag_to_Z_ local_mixing_frac = CS%Gamma_itides, & depth_cutoff = CS%min_zbot_itides*US%Z_to_m) - call read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) + call read_tidal_energy(G, US, tidal_energy_type, tidal_energy_file, CS) !call closeParameterBlock(param_file) @@ -1500,13 +1500,14 @@ end subroutine post_tidal_diagnostics ! TODO: move this subroutine to MOM_internal_tide_input module (?) !> This subroutine read tidal energy inputs from a file. -subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) +subroutine read_tidal_energy(G, US, tidal_energy_type, tidal_energy_file, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type character(len=20), intent(in) :: tidal_energy_type !< The type of tidal energy inputs to read character(len=200), intent(in) :: tidal_energy_file !< The file from which to read tidalinputs type(tidal_mixing_cs), pointer :: CS !< The control structure for this module ! local - integer :: isd, ied, jsd, jed, nz + integer :: i, j, isd, ied, jsd, jed, nz real, allocatable, dimension(:,:) :: tidal_energy_flux_2d ! input tidal energy flux at T-grid points (W/m^2) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke @@ -1516,10 +1517,12 @@ subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed)) allocate(tidal_energy_flux_2d(isd:ied,jsd:jed)) call MOM_read_data(tidal_energy_file,'wave_dissipation',tidal_energy_flux_2d, G%domain) - CS%tidal_qe_2d = (CS%Gamma_itides) * tidal_energy_flux_2d + do j=G%jsc,G%jec ; do i=G%isc,G%iec + CS%tidal_qe_2d(i,j) = CS%Gamma_itides * tidal_energy_flux_2d(i,j) + enddo ; enddo deallocate(tidal_energy_flux_2d) case ('ER03') ! Egbert & Ray 2003 - call read_tidal_constituents(G, tidal_energy_file, CS) + call read_tidal_constituents(G, US, tidal_energy_file, CS) case default call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.") end select @@ -1527,31 +1530,33 @@ subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) end subroutine read_tidal_energy !> This subroutine reads tidal input energy from a file by constituent. -subroutine read_tidal_constituents(G, tidal_energy_file, CS) +subroutine read_tidal_constituents(G, US, tidal_energy_file, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type character(len=200), intent(in) :: tidal_energy_file !< The file from which to read tidal energy inputs type(tidal_mixing_cs), pointer :: CS !< The control structure for this module - ! local - integer :: k, isd, ied, jsd, jed, i,j - integer, dimension(4) :: nz_in - real, parameter :: p33 = 1.0/3.0 + ! local variables + real, parameter :: C1_3 = 1.0/3.0 real, dimension(SZI_(G),SZJ_(G)) :: & tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert real, allocatable, dimension(:) :: & - z_t, & ! depth from surface to midpoint of input layer [cm] - z_w ! depth from surface to top of input layer [cm] + z_t, & ! depth from surface to midpoint of input layer [Z] + z_w ! depth from surface to top of input layer [Z] real, allocatable, dimension(:,:,:) :: & tc_m2, & ! input lunar semidiurnal tidal energy flux [W/m^2] tc_s2, & ! input solar semidiurnal tidal energy flux [W/m^2] tc_k1, & ! input lunar diurnal tidal energy flux [W/m^2] tc_o1 ! input lunar diurnal tidal energy flux [W/m^2] + integer, dimension(4) :: nz_in + integer :: k, is, ie, js, je, isd, ied, jsd, jed, i, j + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ! get number of input levels: - call field_size(tidal_energy_file, 'z_t',nz_in) + call field_size(tidal_energy_file, 'z_t', nz_in) ! allocate local variables allocate(z_t(nz_in(1)), z_w(nz_in(1)) ) @@ -1569,28 +1574,30 @@ subroutine read_tidal_constituents(G, tidal_energy_file, CS) call MOM_read_data(tidal_energy_file, 'S2', tc_s2, G%domain) call MOM_read_data(tidal_energy_file, 'K1', tc_k1, G%domain) call MOM_read_data(tidal_energy_file, 'O1', tc_o1, G%domain) - call MOM_read_data(tidal_energy_file, 'z_t', z_t) - call MOM_read_data(tidal_energy_file, 'z_w', z_w) - - !### THE USE OF WHERE STTAEMENTS IS STRONGLY DISCOURAGED IN MOM6! - where (abs(G%geoLatT(:,:)) < 30.0) - tidal_qk1(:,:) = p33 - tidal_qo1(:,:) = p33 - elsewhere - tidal_qk1(:,:) = 1.0 - tidal_qo1(:,:) = 1.0 - endwhere - - CS%tidal_qe_3d_in = 0.0 + ! Note the hard-coded assumption that z_t and z_w in the file are in centimeters. + call MOM_read_data(tidal_energy_file, 'z_t', z_t, scale=100.0*US%m_to_Z) + call MOM_read_data(tidal_energy_file, 'z_w', z_w, scale=100.0*US%m_to_Z) + + do j=js,je ; do i=is,ie + if (abs(G%geoLatT(i,j)) < 30.0) then + tidal_qk1(i,j) = C1_3 + tidal_qo1(i,j) = C1_3 + else + tidal_qk1(i,j) = 1.0 + tidal_qo1(i,j) = 1.0 + endif + enddo ; enddo + + CS%tidal_qe_3d_in(:,:,:) = 0.0 do k=1,nz_in(1) - ! input cell thickness - CS%h_src(k) = (z_t(k)-z_w(k))*2.0 *1e-2 + ! Store the input cell thickness in m for use with CVmix. + CS%h_src(k) = US%Z_to_m*(z_t(k)-z_w(k))*2.0 ! form tidal_qe_3d_in from weighted tidal constituents - !### THE USE OF WHERE STATEMENTS IS STRONGLY DISCOURAGED IN MOM6! - where (((z_t(k)*1e-2) <= G%bathyT(:,:)*G%Zd_to_m) .and. (z_w(k)*1e-2 > CS%tidal_diss_lim_tc)) - CS%tidal_qe_3d_in(:,:,k) = p33*tc_m2(:,:,k) + p33*tc_s2(:,:,k) + & - tidal_qk1*tc_k1(:,:,k) + tidal_qo1*tc_o1(:,:,k) - endwhere + do j=js,je ; do i=is,ie + if ((z_t(k) <= G%bathyT(i,j)) .and. (z_w(k) > CS%tidal_diss_lim_tc)) & + CS%tidal_qe_3d_in(i,j,k) = C1_3*tc_m2(i,j,k) + C1_3*tc_s2(i,j,k) + & + tidal_qk1(i,j)*tc_k1(i,j,k) + tidal_qo1(i,j)*tc_o1(i,j,k) + enddo ; enddo enddo !open(unit=1905,file="out_1905.txt",access="APPEND") @@ -1601,7 +1608,7 @@ subroutine read_tidal_constituents(G, tidal_energy_file, CS) ! do k=50,nz_in(1) ! write(1905,*) i,j,k ! write(1905,*) CS%tidal_qe_3d_in(i,j,k), tc_m2(i,j,k) - ! write(1905,*) z_t(k), G%bathyT(i,j)*G%Zd_to_m, z_w(k),CS%tidal_diss_lim_tc + ! write(1905,*) z_t(k), G%bathyT(i,j), z_w(k),CS%tidal_diss_lim_tc ! end do ! endif ! enddo @@ -1614,12 +1621,11 @@ subroutine read_tidal_constituents(G, tidal_energy_file, CS) endif !! collapse 3D q*E to 2D q*E - !CS%tidal_qe_2d = 0.0 - !do k=1,nz_in(1) - ! where (z_t(k) <= G%bathyT(:,:)*G%Zd_to_m) - ! CS%tidal_qe_2d(:,:) = CS%tidal_qe_2d(:,:) + CS%tidal_qe_3d_in(:,:,k) - ! endwhere - !enddo + !CS%tidal_qe_2d(:,:) = 0.0 + !do k=1,nz_in(1) ; do j=js,je ; do i=is,ie + ! if (z_t(k) <= G%bathyT(i,j)) & + ! CS%tidal_qe_2d(i,j) = CS%tidal_qe_2d(i,j) + CS%tidal_qe_3d_in(i,j,k) + !enddo ; enddo ; enddo ! initialize input remapping: call initialize_remapping(CS%remap_cs, remapping_scheme="PLM", & From 9ac67ccaa1f64f8bda30c03704bb92adbe4232e5 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 14 Nov 2018 17:53:41 -0500 Subject: [PATCH 04/11] +Add unit_scale_type argument to tracer_Z_init Use US%m_to_Z in place of 1/G%Zd_to_m to convert units from m to Z in tracer_Z_init. This required adding a unit_scale_type argument to tracer_Z_init and 4 subroutines that used the previous interface to tracer_Z_init. All answers are bitwise identical. --- src/tracer/MOM_OCMIP2_CFC.F90 | 15 +++++++++------ src/tracer/MOM_generic_tracer.F90 | 8 +++++--- src/tracer/MOM_tracer_Z_init.F90 | 6 ++++-- src/tracer/MOM_tracer_flow_control.F90 | 8 ++++---- src/tracer/ideal_age_example.F90 | 8 +++++--- src/tracer/oil_tracer.F90 | 11 ++++++----- 6 files changed, 33 insertions(+), 23 deletions(-) diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90 index ebff38508c..6712088988 100644 --- a/src/tracer/MOM_OCMIP2_CFC.F90 +++ b/src/tracer/MOM_OCMIP2_CFC.F90 @@ -18,6 +18,7 @@ module MOM_OCMIP2_CFC use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type @@ -308,13 +309,14 @@ subroutine flux_init_OCMIP2_CFC(CS, verbosity) end subroutine flux_init_OCMIP2_CFC !> Initialize the OCMP2 CFC tracer fields and set up the tracer output. -subroutine initialize_OCMIP2_CFC(restart, day, G, GV, h, diag, OBC, CS, & +subroutine initialize_OCMIP2_CFC(restart, day, G, GV, US, h, diag, OBC, CS, & sponge_CSp, diag_to_Z_CSp) logical, intent(in) :: restart !< .true. if the fields have already been !! read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thicknesses, in H !! (usually m or kg m-2). @@ -343,12 +345,12 @@ subroutine initialize_OCMIP2_CFC(restart, day, G, GV, h, diag, OBC, CS, & if (.not.restart .or. (CS%tracers_may_reinit .and. & .not.query_initialized(CS%CFC11, CS%CFC11_name, CS%restart_CSp))) & call init_tracer_CFC(h, CS%CFC11, CS%CFC11_name, CS%CFC11_land_val, & - CS%CFC11_IC_val, G, CS) + CS%CFC11_IC_val, G, US, CS) if (.not.restart .or. (CS%tracers_may_reinit .and. & .not.query_initialized(CS%CFC12, CS%CFC12_name, CS%restart_CSp))) & call init_tracer_CFC(h, CS%CFC12, CS%CFC12_name, CS%CFC12_land_val, & - CS%CFC12_IC_val, G, CS) + CS%CFC12_IC_val, G, US, CS) if (associated(OBC)) then ! Steal from updated DOME in the fullness of time. @@ -357,8 +359,9 @@ subroutine initialize_OCMIP2_CFC(restart, day, G, GV, h, diag, OBC, CS, & end subroutine initialize_OCMIP2_CFC !>This subroutine initializes a tracer array. -subroutine init_tracer_CFC(h, tr, name, land_val, IC_val, G, CS) +subroutine init_tracer_CFC(h, tr, name, land_val, IC_val, G, US, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: tr !< The tracer concentration array character(len=*), intent(in) :: name !< The tracer name @@ -378,9 +381,9 @@ subroutine init_tracer_CFC(h, tr, name, land_val, IC_val, G, CS) if (.not.file_exists(CS%IC_file, G%Domain)) & call MOM_error(FATAL, "initialize_OCMIP2_CFC: Unable to open "//CS%IC_file) if (CS%Z_IC_file) then - OK = tracer_Z_init(tr, h, CS%IC_file, name, G) + OK = tracer_Z_init(tr, h, CS%IC_file, name, G, US) if (.not.OK) then - OK = tracer_Z_init(tr, h, CS%IC_file, trim(name), G) + OK = tracer_Z_init(tr, h, CS%IC_file, trim(name), G, US) if (.not.OK) call MOM_error(FATAL,"initialize_OCMIP2_CFC: "//& "Unable to read "//trim(name)//" from "//& trim(CS%IC_file)//".") diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index ee1f038180..12f7ecf0c4 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -40,6 +40,7 @@ module MOM_generic_tracer use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_Z_init, only : tracer_Z_init use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z + use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs use MOM_open_boundary, only : ocean_OBC_type use MOM_verticalGrid, only : verticalGrid_type @@ -222,13 +223,14 @@ end function register_MOM_generic_tracer !! !! This subroutine initializes the NTR tracer fields in tr(:,:,:,:) !! and it sets up the tracer output. - subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, diag, OBC, CS, & + subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, CS, & sponge_CSp, ALE_sponge_CSp,diag_to_Z_CSp) logical, intent(in) :: restart !< .true. if the fields have already been !! read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(diag_ctrl), target, intent(in) :: diag !< Regulates diagnostic output. @@ -319,9 +321,9 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, & "initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file) if (CS%Z_IC_file) then - OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G) + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G, US) if (.not.OK) then - OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G) + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G, US) if (.not.OK) call MOM_error(FATAL,"initialize_MOM_Generic_tracer: "//& "Unable to read "//trim(g_tracer_name)//" from "//& trim(CS%IC_file)//".") diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 7450571500..d78821cd46 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -8,6 +8,7 @@ module MOM_tracer_Z_init ! use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_io, only : MOM_read_data +use MOM_unit_scaling, only : unit_scale_type use netcdf @@ -21,9 +22,10 @@ module MOM_tracer_Z_init !> This function initializes a tracer by reading a Z-space file, returning !! .true. if this appears to have been successful, and false otherwise. -function tracer_Z_init(tr, h, filename, tr_name, G, missing_val, land_val) +function tracer_Z_init(tr, h, filename, tr_name, G, US, missing_val, land_val) logical :: tracer_Z_init !< A return code indicating if the initialization has been successful type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(out) :: tr !< The tracer to initialize real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & @@ -82,7 +84,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, missing_val, land_val) ! Find out the number of input levels and read the depth of the edges, ! also modifying their sign convention to be monotonically decreasing. call read_Z_edges(filename, tr_name, z_edges, nz_in, has_edges, use_missing, & - missing, scale=1.0/G%Zd_to_m) + missing, scale=US%m_to_Z) if (nz_in < 1) then tracer_Z_init = .false. return diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 3dcd9ac192..6438a55ed2 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -317,23 +317,23 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag call initialize_ISOMIP_tracer(restart, day, G, GV, h, diag, OBC, CS%ISOMIP_tracer_CSp, & ALE_sponge_CSp, diag_to_Z_CSp) if (CS%use_ideal_age) & - call initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS%ideal_age_tracer_CSp, & + call initialize_ideal_age_tracer(restart, day, G, GV, US, h, diag, OBC, CS%ideal_age_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) if (CS%use_regional_dyes) & call initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS%dye_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) if (CS%use_oil) & - call initialize_oil_tracer(restart, day, G, GV, h, diag, OBC, CS%oil_tracer_CSp, & + call initialize_oil_tracer(restart, day, G, GV, US, h, diag, OBC, CS%oil_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) if (CS%use_advection_test_tracer) & call initialize_advection_test_tracer(restart, day, G, GV, h, diag, OBC, CS%advection_test_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) if (CS%use_OCMIP2_CFC) & - call initialize_OCMIP2_CFC(restart, day, G, GV, h, diag, OBC, CS%OCMIP2_CFC_CSp, & + call initialize_OCMIP2_CFC(restart, day, G, GV, US, h, diag, OBC, CS%OCMIP2_CFC_CSp, & sponge_CSp, diag_to_Z_CSp) #ifdef _USE_GENERIC_TRACER if (CS%use_MOM_generic_tracer) & - call initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, diag, OBC, & + call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, & CS%MOM_generic_tracer_CSp, sponge_CSp, ALE_sponge_CSp, diag_to_Z_CSp) #endif if (CS%use_pseudo_salt_tracer) & diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index d7fcb53324..750fa83021 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/src/tracer/ideal_age_example.F90 @@ -18,6 +18,7 @@ module ideal_age_example use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type @@ -193,13 +194,14 @@ function register_ideal_age_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) end function register_ideal_age_tracer !> Sets the ideal age traces to their initial values and sets up the tracer output -subroutine initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS, & +subroutine initialize_ideal_age_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & sponge_CSp, diag_to_Z_CSp) logical, intent(in) :: restart !< .true. if the fields have already !! been read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate @@ -250,10 +252,10 @@ subroutine initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS, & if (CS%Z_IC_file) then OK = tracer_Z_init(CS%tr(:,:,:,m), h, CS%IC_file, name,& - G, -1e34, 0.0) ! CS%land_val(m)) + G, US, -1e34, 0.0) ! CS%land_val(m)) if (.not.OK) then OK = tracer_Z_init(CS%tr(:,:,:,m), h, CS%IC_file, & - trim(name), G, -1e34, 0.0) ! CS%land_val(m)) + trim(name), G, US, -1e34, 0.0) ! CS%land_val(m)) if (.not.OK) call MOM_error(FATAL,"initialize_ideal_age_tracer: "//& "Unable to read "//trim(name)//" from "//& trim(CS%IC_file)//".") diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90 index 3b98c19a73..3130ba3804 100644 --- a/src/tracer/oil_tracer.F90 +++ b/src/tracer/oil_tracer.F90 @@ -18,8 +18,8 @@ module oil_tracer use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init -use MOM_variables, only : surface -use MOM_variables, only : thermo_var_ptrs +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use coupler_types_mod, only : coupler_type_set_data, ind_csurf @@ -201,13 +201,14 @@ function register_oil_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) end function register_oil_tracer !> Initialize the oil tracers and set up tracer output -subroutine initialize_oil_tracer(restart, day, G, GV, h, diag, OBC, CS, & +subroutine initialize_oil_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & sponge_CSp, diag_to_Z_CSp) logical, intent(in) :: restart !< .true. if the fields have already !! been read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate @@ -266,10 +267,10 @@ subroutine initialize_oil_tracer(restart, day, G, GV, h, diag, OBC, CS, & if (CS%Z_IC_file) then OK = tracer_Z_init(CS%tr(:,:,:,m), h, CS%IC_file, name, & - G, -1e34, 0.0) ! CS%land_val(m)) + G, US, -1e34, 0.0) ! CS%land_val(m)) if (.not.OK) then OK = tracer_Z_init(CS%tr(:,:,:,m), h, CS%IC_file, & - trim(name), G, -1e34, 0.0) ! CS%land_val(m)) + trim(name), G, US, -1e34, 0.0) ! CS%land_val(m)) if (.not.OK) call MOM_error(FATAL,"initialize_oil_tracer: "//& "Unable to read "//trim(name)//" from "//& trim(CS%IC_file)//".") From d2bfd736dfc6f734b653dbae1550da8979b20a61 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 14 Nov 2018 18:58:19 -0500 Subject: [PATCH 05/11] +Add unit_scale_type argument to diag_remap_update Added a unit_scale_type argument to diag_remap_update and a pointer to a unit_scale_type structure to diag_mediator_init and store a this pointer in the diag_ctrl type, all to accomodate rescaling of depths via US%Z_to_m instead of G%Zd_to_m. All answers are bitwise identical. --- src/core/MOM.F90 | 2 +- src/framework/MOM_diag_mediator.F90 | 9 ++++++--- src/framework/MOM_diag_remap.F90 | 9 +++++---- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9bb2fa74cf..3126b0b8ec 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2256,7 +2256,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & diag => CS%diag ! Initialize the diag mediator. - call diag_mediator_init(G, GV, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory) + call diag_mediator_init(G, GV, US, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory) if (present(diag_ptr)) diag_ptr => CS%diag ! Initialize the diagnostics masks for native arrays. diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 117adb512e..db9e391610 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -15,7 +15,7 @@ module MOM_diag_mediator use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type -use MOM_unit_scaling, only : unit_scale_type +use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : EOS_type use MOM_diag_remap, only : diag_remap_ctrl @@ -223,6 +223,7 @@ module MOM_diag_mediator type(EOS_type), pointer :: eqn_of_state => null() !< The equation of state type type(ocean_grid_type), pointer :: G => null() !< The ocean grid type type(verticalGrid_type), pointer :: GV => null() !< The model's vertical ocean grid + type(unit_scale_type), pointer :: US => null() !< A dimensional unit scaling type !> The volume cell measure (special diagnostic) manager id integer :: volume_cell_measure_dm_id = -1 @@ -2177,9 +2178,10 @@ end subroutine diag_mediator_infrastructure_init !> diag_mediator_init initializes the MOM diag_mediator and opens the available !! diagnostics file, if appropriate. -subroutine diag_mediator_init(G, GV, nz, param_file, diag_cs, doc_file_dir) +subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) type(ocean_grid_type), target, intent(inout) :: G !< The ocean grid type. type(verticalGrid_type), target, intent(in) :: GV !< The ocean vertical grid structure + type(unit_scale_type), target, intent(in) :: US !< A dimensional unit scaling type integer, intent(in) :: nz !< The number of layers in the model's native grid. type(param_file_type), intent(in) :: param_file !< Parameter file structure type(diag_ctrl), intent(inout) :: diag_cs !< A pointer to a type with many variables @@ -2251,6 +2253,7 @@ subroutine diag_mediator_init(G, GV, nz, param_file, diag_cs, doc_file_dir) ! Keep pointers grid, h, T, S needed diagnostic remapping diag_cs%G => G diag_cs%GV => GV + diag_cs%US => US diag_cs%h => null() diag_cs%T => null() diag_cs%S => null() @@ -2404,7 +2407,7 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S) do i=1, diag_cs%num_diag_coords call diag_remap_update(diag_cs%diag_remap_cs(i), & - diag_cs%G, diag_cs%GV, h_diag, T_diag, S_diag, & + diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & diag_cs%eqn_of_state) enddo diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 748c09ba0c..632258d5d2 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -225,10 +225,11 @@ function diag_remap_axes_configured(remap_cs) !! height or layer thicknesses changes. In the case of density-based !! coordinates then technically we should also regenerate the !! target grid whenever T/S change. -subroutine diag_remap_update(remap_cs, G, GV, h, T, S, eqn_of_state) +subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure type(ocean_grid_type), pointer :: G !< The ocean's grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(:, :, :), intent(in) :: h !< New thickness real, dimension(:, :, :), intent(in) :: T !< New T real, dimension(:, :, :), intent(in) :: S !< New S @@ -278,15 +279,15 @@ subroutine diag_remap_update(remap_cs, G, GV, h, T, S, eqn_of_state) GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then call build_rho_column(get_rho_CS(remap_cs%regrid_cs), G%ke, & - G%Zd_to_m*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & + US%Z_to_m*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect, h_neglect_edge) elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then ! call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, & -! G%Zd_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) +! US%Z_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) call MOM_error(FATAL,"diag_remap_update: SLIGHT coordinate not coded for diagnostics yet!") elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then ! call build_hycom1_column(remap_cs%regrid_cs, nz, & -! G%Zd_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) +! US%Z_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") endif remap_cs%h(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1) From 366efc8c6e6ddb7c337d0ab565291d68e81dc9c0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Nov 2018 04:45:20 -0500 Subject: [PATCH 06/11] +Add unit_scale_type arg to MOM_initialize_fixed Added a unit_scale_type argument to MOM_initialize_fixed and moved the call to rescale_dyn_horgrid_bathymetry into MOM_initialize_fixed immediately after the call to MOM_initialize_topography. Also added a unit_scale_type argument to mask_outside_OBCs and open_boundary_config to accomodate rescaling of depths via US%Z_to_m instead of G%Zd_to_m. All answers are bitwise identical. --- src/core/MOM.F90 | 7 ++---- src/core/MOM_open_boundary.F90 | 24 ++++++++++-------- .../MOM_fixed_initialization.F90 | 25 +++++++++++-------- 3 files changed, 29 insertions(+), 27 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3126b0b8ec..4063ac619d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -72,7 +72,6 @@ module MOM use MOM_dynamics_unsplit_RK2, only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2 use MOM_dynamics_unsplit_RK2, only : MOM_dyn_unsplit_RK2_CS use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid -use MOM_dyn_horgrid, only : rescale_dyn_horgrid_bathymetry use MOM_EOS, only : EOS_init, calculate_density, calculate_TFreeze use MOM_fixed_initialization, only : MOM_initialize_fixed use MOM_grid, only : ocean_grid_type, MOM_grid_init, MOM_grid_end @@ -1978,11 +1977,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call MOM_timing_init(CS) ! Allocate initialize time-invariant MOM variables. - call MOM_initialize_fixed(dG, CS%OBC, param_file, write_geom_files, dirs%output_directory) + call MOM_initialize_fixed(dG, US, CS%OBC, param_file, write_geom_files, dirs%output_directory) call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)") - if (dG%Zd_to_m /= US%Z_to_m) call rescale_dyn_horgrid_bathymetry(dG, US%Z_to_m) - if (associated(CS%OBC)) call call_OBC_register(param_file, CS%update_OBC_CSp, CS%OBC) call tracer_registry_init(param_file, CS%tracer_Reg) @@ -2982,7 +2979,7 @@ subroutine extract_surface_state(CS, sfc_state) numberOfErrors=0 ! count number of errors do j=js,je; do i=is,ie if (G%mask2dT(i,j)>0.) then - bathy_m = G%Zd_to_m*G%bathyT(i,j) + bathy_m = CS%US%Z_to_m * G%bathyT(i,j) localError = sfc_state%sea_lev(i,j)<=-bathy_m & .or. sfc_state%sea_lev(i,j)>= CS%bad_val_ssh_max & .or. sfc_state%sea_lev(i,j)<=-CS%bad_val_ssh_max & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 6296dbc35b..8fd4cc06f5 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -19,12 +19,13 @@ module MOM_open_boundary use MOM_obsolete_params, only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char use MOM_string_functions, only : extract_word, remove_spaces use MOM_tracer_registry, only : tracer_type, tracer_registry_type, tracer_name_lookup -use MOM_variables, only : thermo_var_ptrs use time_interp_external_mod, only : init_external_field, time_interp_external use time_interp_external_mod, only : time_interp_external_init use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping use MOM_regridding, only : regridding_CS +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -290,8 +291,9 @@ module MOM_open_boundary !> here. The remainder of the segment data are initialized in a !> later call to update_open_boundary_data -subroutine open_boundary_config(G, param_file, OBC) +subroutine open_boundary_config(G, US, param_file, OBC) type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file handle type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure ! Local variables @@ -473,7 +475,7 @@ subroutine open_boundary_config(G, param_file, OBC) "is entering the domain.", units="m", default=0.0) endif - if (mask_outside) call mask_outside_OBCs(G, param_file, OBC) + if (mask_outside) call mask_outside_OBCs(G, US, param_file, OBC) ! All tracers are using the same restoring length scale for now, but we may want to make this ! tracer-specific in the future for example, in cases where certain tracers are poorly constrained @@ -3690,17 +3692,18 @@ end subroutine fill_temp_salt_segments !> Find the region outside of all open boundary segments and !! make sure it is set to land mask. Gonna need to know global land !! mask as well to get it right... -subroutine mask_outside_OBCs(G, param_file, OBC) - type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure - type(param_file_type), intent(in) :: param_file !< Parameter file handle - type(ocean_OBC_type), pointer :: OBC !< Open boundary structure +subroutine mask_outside_OBCs(G, US, param_file, OBC) + type(dyn_horgrid_type), intent(inout) :: G !< Ocean grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file handle + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type -! Local variables + ! Local variables integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n integer :: i, j logical :: fatal_error = .False. real :: min_depth - integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 + integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 character(len=256) :: mesg ! Message for error messages. type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list real, allocatable, dimension(:,:) :: color, color2 ! For sorting inside from outside, @@ -3709,8 +3712,7 @@ subroutine mask_outside_OBCs(G, param_file, OBC) if (.not. associated(OBC)) return call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - default=0.0, do_not_log=.true.) - min_depth = min_depth / G%Zd_to_m + units="m", default=0.0, scale=US%m_to_Z, do_not_log=.true.) allocate(color(G%isd:G%ied, G%jsd:G%jed)) ; color = 0 allocate(color2(G%isd:G%ied, G%jsd:G%jed)) ; color2 = 0 diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index b754b19bcb..24d496703c 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -6,7 +6,7 @@ module MOM_fixed_initialization use MOM_debugging, only : hchksum, qchksum, uvchksum use MOM_domains, only : pass_var -use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_dyn_horgrid, only : dyn_horgrid_type, rescale_dyn_horgrid_bathymetry use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type @@ -25,6 +25,7 @@ module MOM_fixed_initialization use MOM_shared_initialization, only : reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list use MOM_shared_initialization, only : read_face_length_list, set_velocity_depth_max, set_velocity_depth_min use MOM_shared_initialization, only : compute_global_grid_integrals, write_ocean_geometry_file +use MOM_unit_scaling, only : unit_scale_type use user_initialization, only : user_initialize_topography use DOME_initialization, only : DOME_initialize_topography @@ -51,8 +52,9 @@ module MOM_fixed_initialization ! ----------------------------------------------------------------------------- !> MOM_initialize_fixed sets up time-invariant quantities related to MOM6's !! horizontal grid, bathymetry, and the Coriolis parameter. -subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) +subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) type(dyn_horgrid_type), intent(inout) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ocean_OBC_type), pointer :: OBC !< Open boundary structure. type(param_file_type), intent(in) :: PF !< A structure indicating the open file !! to parse for model parameter values. @@ -75,19 +77,20 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) "The directory in which input files are found.", default=".") inputdir = slasher(inputdir) -! Set up the parameters of the physical domain (i.e. the grid), G + ! Set up the parameters of the physical domain (i.e. the grid), G call set_grid_metrics(G, PF) -! Set up the bottom depth, G%bathyT either analytically or from file -! This also sets G%max_depth based on the input parameter MAXIMUM_DEPTH, -! or, if absent, is diagnosed as G%max_depth = max( G%D(:,:) ) + ! Set up the bottom depth, G%bathyT either analytically or from file + ! This also sets G%max_depth based on the input parameter MAXIMUM_DEPTH, + ! or, if absent, is diagnosed as G%max_depth = max( G%D(:,:) ) call MOM_initialize_topography(G%bathyT, G%max_depth, G, PF) + if (G%Zd_to_m /= US%Z_to_m) call rescale_dyn_horgrid_bathymetry(G, US%Z_to_m) ! To initialize masks, the bathymetry in halo regions must be filled in call pass_var(G%bathyT, G%Domain) -! Determine the position of any open boundaries - call open_boundary_config(G, PF, OBC) + ! Determine the position of any open boundaries + call open_boundary_config(G, US, PF, OBC) ! Make bathymetry consistent with open boundaries call open_boundary_impose_normal_slope(OBC, G, G%bathyT) @@ -99,14 +102,14 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) call open_boundary_impose_land_mask(OBC, G, G%areaCu, G%areaCv) if (debug) then - call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1) + call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1, scale=US%Z_to_m) call hchksum(G%mask2dT, 'MOM_initialize_fixed: mask2dT ', G%HI) call uvchksum('MOM_initialize_fixed: mask2dC[uv]', G%mask2dCu, & G%mask2dCv, G%HI) call qchksum(G%mask2dBu, 'MOM_initialize_fixed: mask2dBu ', G%HI) endif -! Modulate geometric scales according to geography. + ! Modulate geometric scales according to geography. call get_param(PF, mdl, "CHANNEL_CONFIG", config, & "A parameter that determines which set of channels are \n"//& "restricted to specific widths. Options are:\n"//& @@ -128,7 +131,7 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) "Unrecognized channel configuration "//trim(config)) end select -! This call sets the topography at velocity points. + ! This call sets the topography at velocity points. if (G%bathymetry_at_vel) then call get_param(PF, mdl, "VELOCITY_DEPTH_CONFIG", config, & "A string that determines how the topography is set at \n"//& From 6c478014368a4e8845a6c71f6c086cd011451efc Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Nov 2018 06:02:02 -0500 Subject: [PATCH 07/11] +Add unit_scale_type arg to MOM_sum_output_init Added a unit_scale_type argument to MOM_sum_output_init, depth_list_setup, read_depth_list and write_depth_list, and used elements of this type to rescale depth variables in place of G%Zd_to_m. All answers are bitwise identical. --- src/core/MOM.F90 | 2 +- src/diagnostics/MOM_sum_output.F90 | 32 ++++++++++++++++-------------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4063ac619d..b2d211796f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2452,7 +2452,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%nstep_tot = 0 if (present(count_calls)) CS%count_calls = count_calls - call MOM_sum_output_init(G, param_file, dirs%output_directory, & + call MOM_sum_output_init(G, US, param_file, dirs%output_directory, & CS%ntrunc, Time_init, CS%sum_output_CSp) ! Flag whether to save initial conditions in finish_MOM_initialization() or not. diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 112c3b9104..92028dabf4 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -117,9 +117,10 @@ module MOM_sum_output contains !> MOM_sum_output_init initializes the parameters and settings for the MOM_sum_output module. -subroutine MOM_sum_output_init(G, param_file, directory, ntrnc, & +subroutine MOM_sum_output_init(G, US, param_file, directory, ntrnc, & Input_start_time, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time !! parameters. character(len=*), intent(in) :: directory !< The directory where the energy file goes. @@ -213,7 +214,7 @@ subroutine MOM_sum_output_init(G, param_file, directory, ntrnc, & call get_param(param_file, mdl, "DEPTH_LIST_MIN_INC", CS%D_list_min_inc, & "The minimum increment between the depths of the \n"//& "entries in the depth-list file.", & - units="m", default=1.0E-10, scale=1.0/G%Zd_to_m) + units="m", default=1.0E-10, scale=US%m_to_Z) if (CS%read_depth_list) then call get_param(param_file, mdl, "DEPTH_LIST_FILE", CS%depth_list_file, & "The name of the depth list file.", default="Depth_list.nc") @@ -222,7 +223,7 @@ subroutine MOM_sum_output_init(G, param_file, directory, ntrnc, & endif allocate(CS%lH(G%ke)) - call depth_list_setup(G, CS) + call depth_list_setup(G, US, CS) else CS%list_size = 0 endif @@ -1023,8 +1024,9 @@ end subroutine accumulate_net_input !! cross sectional areas at each depth and the volume of fluid deeper !! than each depth. This might be read from a previously created file !! or it might be created anew. (For now only new creation occurs. -subroutine depth_list_setup(G, CS) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure +subroutine depth_list_setup(G, US, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(Sum_output_CS), pointer :: CS !< The control structure returned by a !! previous call to MOM_sum_output_init. ! Local variables @@ -1032,13 +1034,13 @@ subroutine depth_list_setup(G, CS) if (CS%read_depth_list) then if (file_exists(CS%depth_list_file)) then - call read_depth_list(G, CS, CS%depth_list_file) + call read_depth_list(G, US, CS, CS%depth_list_file) else if (is_root_pe()) call MOM_error(WARNING, "depth_list_setup: "// & trim(CS%depth_list_file)//" does not exist. Creating a new file.") call create_depth_list(G, CS) - call write_depth_list(G, CS, CS%depth_list_file, CS%list_size+1) + call write_depth_list(G, US, CS, CS%depth_list_file, CS%list_size+1) endif else call create_depth_list(G, CS) @@ -1177,8 +1179,9 @@ subroutine create_depth_list(G, CS) end subroutine create_depth_list !> This subroutine writes out the depth list to the specified file. -subroutine write_depth_list(G, CS, filename, list_size) +subroutine write_depth_list(G, US, CS, filename, list_size) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(Sum_output_CS), pointer :: CS !< The control structure returned by a !! previous call to MOM_sum_output_init. character(len=*), intent(in) :: filename !< The path to the depth list file to write. @@ -1235,7 +1238,7 @@ subroutine write_depth_list(G, CS, filename, list_size) if (status /= NF90_NOERR) call MOM_error(WARNING, & filename//trim(NF90_STRERROR(status))) - do k=1,list_size ; tmp(k) = G%Zd_to_m*CS%DL(k)%depth ; enddo + do k=1,list_size ; tmp(k) = US%Z_to_m*CS%DL(k)%depth ; enddo status = NF90_PUT_VAR(ncid, Did, tmp) if (status /= NF90_NOERR) call MOM_error(WARNING, & filename//" depth "//trim(NF90_STRERROR(status))) @@ -1245,7 +1248,7 @@ subroutine write_depth_list(G, CS, filename, list_size) if (status /= NF90_NOERR) call MOM_error(WARNING, & filename//" area "//trim(NF90_STRERROR(status))) - do k=1,list_size ; tmp(k) = G%Zd_to_m*CS%DL(k)%vol_below ; enddo + do k=1,list_size ; tmp(k) = US%Z_to_m*CS%DL(k)%vol_below ; enddo status = NF90_PUT_VAR(ncid, Vid, tmp) if (status /= NF90_NOERR) call MOM_error(WARNING, & filename//" vol_below "//trim(NF90_STRERROR(status))) @@ -1258,8 +1261,9 @@ end subroutine write_depth_list !> This subroutine reads in the depth list to the specified file !! and allocates and sets up CS%DL and CS%list_size . -subroutine read_depth_list(G, CS, filename) +subroutine read_depth_list(G, US, CS, filename) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(Sum_output_CS), pointer :: CS !< The control structure returned by a !! previous call to MOM_sum_output_init. character(len=*), intent(in) :: filename !< The path to the depth list file to read. @@ -1267,12 +1271,10 @@ subroutine read_depth_list(G, CS, filename) character(len=32) :: mdl character(len=240) :: var_name, var_msg real, allocatable :: tmp(:) - real :: m_to_Z integer :: ncid, status, varid, list_size, k integer :: ndim, len, var_dim_ids(NF90_MAX_VAR_DIMS) mdl = "MOM_sum_output read_depth_list:" - m_to_Z = 1.0/G%Zd_to_m status = NF90_OPEN(filename, NF90_NOWRITE, ncid) if (status /= NF90_NOERR) then @@ -1311,7 +1313,7 @@ subroutine read_depth_list(G, CS, filename) " Difficulties reading variable "//trim(var_msg)//& trim(NF90_STRERROR(status))) - do k=1,list_size ; CS%DL(k)%depth = m_to_Z*tmp(k) ; enddo + do k=1,list_size ; CS%DL(k)%depth = US%m_to_Z*tmp(k) ; enddo var_name = "area" var_msg = trim(var_name)//" in "//trim(filename)//" - " @@ -1337,7 +1339,7 @@ subroutine read_depth_list(G, CS, filename) " Difficulties reading variable "//trim(var_msg)//& trim(NF90_STRERROR(status))) - do k=1,list_size ; CS%DL(k)%vol_below = m_to_Z*tmp(k) ; enddo + do k=1,list_size ; CS%DL(k)%vol_below = US%m_to_Z*tmp(k) ; enddo status = NF90_CLOSE(ncid) if (status /= NF90_NOERR) call MOM_error(WARNING, mdl// & From 3d4c8913de017b510052b389ca3ecbbab4c8bb49 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Nov 2018 06:03:33 -0500 Subject: [PATCH 08/11] +Add optional unit_scale_type arg to initialize_masks Added an optional unit_scale_type argument to initialize_masks and write_ocean_geometry_file, and used elements of this type to rescale depth variables in place of G%Zd_to_m. All answers are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- .../MOM_fixed_initialization.F90 | 6 +++--- src/initialization/MOM_grid_initialize.F90 | 19 ++++++++++++------- .../MOM_shared_initialization.F90 | 19 ++++++++++++------- 4 files changed, 28 insertions(+), 18 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index a7b89c4934..82ae988f4f 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1385,7 +1385,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl ! Set up the bottom depth, G%D either analytically or from file call MOM_initialize_topography(dG%bathyT, G%max_depth, dG, param_file) - if (dG%Zd_to_m /= US%Z_to_m) call rescale_dyn_horgrid_bathymetry(dG, US%Z_to_m) + call rescale_dyn_horgrid_bathymetry(dG, US%Z_to_m) ! Set up the Coriolis parameter, G%f, usually analytically. call MOM_initialize_rotation(dG%CoriolisBu, dG, param_file) diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 24d496703c..c7ebb6a251 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -84,7 +84,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) ! This also sets G%max_depth based on the input parameter MAXIMUM_DEPTH, ! or, if absent, is diagnosed as G%max_depth = max( G%D(:,:) ) call MOM_initialize_topography(G%bathyT, G%max_depth, G, PF) - if (G%Zd_to_m /= US%Z_to_m) call rescale_dyn_horgrid_bathymetry(G, US%Z_to_m) + call rescale_dyn_horgrid_bathymetry(G, US%Z_to_m) ! To initialize masks, the bathymetry in halo regions must be filled in call pass_var(G%bathyT, G%Domain) @@ -96,7 +96,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) call open_boundary_impose_normal_slope(OBC, G, G%bathyT) ! This call sets masks that prohibit flow over any point interpreted as land - call initialize_masks(G, PF) + call initialize_masks(G, PF, US) ! Make OBC mask consistent with land mask call open_boundary_impose_land_mask(OBC, G, G%areaCu, G%areaCv) @@ -162,7 +162,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) call compute_global_grid_integrals(G) ! Write out all of the grid data used by this run. - if (write_geom) call write_ocean_geometry_file(G, PF, output_dir) + if (write_geom) call write_ocean_geometry_file(G, PF, output_dir, US=US) call callTree_leave('MOM_initialize_fixed()') diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 0f5a8505ab..7019564586 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -16,6 +16,7 @@ module MOM_grid_initialize use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_io, only : MOM_read_data, read_data, slasher, file_exists use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE +use MOM_unit_scaling, only : unit_scale_type use mpp_domains_mod, only : mpp_get_domain_extents, mpp_deallocate_domain @@ -1214,26 +1215,30 @@ end function Adcroft_reciprocal !! are 0.0 at any points adjacent to a land point. mask2dBu is 0.0 at !! any land or boundary point. For points in the interior, mask2dCu, !! mask2dCv, and mask2dBu are all 1.0. -subroutine initialize_masks(G, PF) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: PF !< Parameter file structure +subroutine initialize_masks(G, PF, US) + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: PF !< Parameter file structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: Dmin ! The depth for masking in the same units as G%bathyT (Z). - real :: min_depth, mask_depth ! Depths in the same units as G%bathyT (Z). + real :: m_to_Z_scale ! A unit conversion factor from m to Z. + real :: Dmin ! The depth for masking in the same units as G%bathyT (Z). + real :: min_depth ! The minimum ocean depth in the same units as G%bathyT (Z). + real :: mask_depth ! The depth shallower than which to mask a point as land, in Z. character(len=40) :: mdl = "MOM_grid_init initialize_masks" integer :: i, j call callTree_enter("initialize_masks(), MOM_grid_initialize.F90") + m_to_Z_scale = 1.0 ; if (present(US)) m_to_Z_scale = US%m_to_Z call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & "If MASKING_DEPTH is unspecified, then anything shallower than\n"//& "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out.\n"//& "If MASKING_DEPTH is specified, then all depths shallower than\n"//& "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", & - units="m", default=0.0, scale=1.0/G%Zd_to_m) + units="m", default=0.0, scale=m_to_Z_scale) call get_param(PF, mdl, "MASKING_DEPTH", mask_depth, & "The depth below which to mask points as land points, for which all\n"//& "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", & - units="m", default=-9999.0, scale=1.0/G%Zd_to_m) + units="m", default=-9999.0, scale=m_to_Z_scale) Dmin = min_depth if (mask_depth>=0.) Dmin = mask_depth diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 31dfc551b5..9ca352c3ed 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -15,6 +15,7 @@ module MOM_shared_initialization use MOM_io, only : MOM_read_data, MOM_read_vector, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, write_field, var_desc use MOM_string_functions, only : uppercase +use MOM_unit_scaling, only : unit_scale_type use netcdf @@ -1143,12 +1144,13 @@ end subroutine compute_global_grid_integrals ! ----------------------------------------------------------------------------- !> Write out a file describing the topography, Coriolis parameter, grid locations !! and various other fixed fields from the grid. -subroutine write_ocean_geometry_file(G, param_file, directory, geom_file) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid - type(param_file_type), intent(in) :: param_file !< Parameter file structure - character(len=*), intent(in) :: directory !< The directory into which to place the geometry file. - character(len=*), optional, intent(in) :: geom_file !< If present, the name of the geometry file - !! (otherwise the file is "ocean_geometry") +subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid + type(param_file_type), intent(in) :: param_file !< Parameter file structure + character(len=*), intent(in) :: directory !< The directory into which to place the geometry file. + character(len=*), optional, intent(in) :: geom_file !< If present, the name of the geometry file + !! (otherwise the file is "ocean_geometry") + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables. character(len=240) :: filepath @@ -1156,6 +1158,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file) integer, parameter :: nFlds=23 type(vardesc) :: vars(nFlds) type(fieldtype) :: fields(nFlds) + real :: Z_to_m_scale ! A unit conversion factor from Z to m. integer :: unit integer :: file_threading integer :: nFlds_used @@ -1172,6 +1175,8 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + Z_to_m_scale = 1.0 ; if (present(US)) Z_to_m_scale = US%Z_to_m + ! vardesc is a structure defined in MOM_io.F90. The elements of ! this structure, in order, are: ! (1) the variable name for the NetCDF file @@ -1240,7 +1245,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file) call write_field(unit, fields(3), G%Domain%mpp_domain, G%geoLatT) call write_field(unit, fields(4), G%Domain%mpp_domain, G%geoLonT) - do j=js,je ; do i=is,ie ; out_h(i,j) = G%Zd_to_m*G%bathyT(i,j) ; enddo ; enddo + do j=js,je ; do i=is,ie ; out_h(i,j) = Z_to_m_scale*G%bathyT(i,j) ; enddo ; enddo call write_field(unit, fields(5), G%Domain%mpp_domain, out_h) call write_field(unit, fields(6), G%Domain%mpp_domain, G%CoriolisBu) From fac1464e831ebe2b6165e27dc2fc014023fee1bd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Nov 2018 10:23:13 -0500 Subject: [PATCH 09/11] +Eliminated Zd_to_m from grid types Eliminated Zd_to_m from dyn_horgrid_type and ocean_grid_type. Instead, any dimensional rescaling uses elements of unit_scale_types. All answers are bitwise identical. --- src/core/MOM_grid.F90 | 8 +++----- src/core/MOM_transcribe_grid.F90 | 2 -- src/framework/MOM_dyn_horgrid.F90 | 16 +++++++--------- 3 files changed, 10 insertions(+), 16 deletions(-) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index c0ca264d68..e226598b7f 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -132,7 +132,6 @@ module MOM_grid real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & bathyT !< Ocean bottom depth at tracer points, in depth units. - real :: Zd_to_m = 1.0 !< The conversion factor between the units of bathyT and m. logical :: bathymetry_at_vel !< If true, there are separate values for the !! basin depths at velocity points. Otherwise the effects of @@ -165,7 +164,7 @@ module MOM_grid real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain real :: Rad_Earth = 6.378e6 !< The radius of the planet in meters. - real :: max_depth !< The maximum depth of the ocean in depth units (scaled by Zd_to_m). + real :: max_depth !< The maximum depth of the ocean in depth units (Z). end type ocean_grid_type contains @@ -359,13 +358,13 @@ subroutine rescale_grid_bathymetry(G, m_in_new_units) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - if (m_in_new_units == G%Zd_to_m) return + if (m_in_new_units == 1.0) return if (m_in_new_units < 0.0) & call MOM_error(FATAL, "rescale_grid_bathymetry: Negative depth units are not permitted.") if (m_in_new_units == 0.0) & call MOM_error(FATAL, "rescale_grid_bathymetry: Zero depth units are not permitted.") - rescale = G%Zd_to_m / m_in_new_units + rescale = 1.0 / m_in_new_units do j=jsd,jed ; do i=isd,ied G%bathyT(i,j) = rescale*G%bathyT(i,j) enddo ; enddo @@ -376,7 +375,6 @@ subroutine rescale_grid_bathymetry(G, m_in_new_units) G%Dblock_v(i,J) = rescale*G%Dblock_v(i,J) ; G%Dopen_v(i,J) = rescale*G%Dopen_v(i,J) enddo ; enddo ; endif G%max_depth = rescale*G%max_depth - G%Zd_to_m = m_in_new_units end subroutine rescale_grid_bathymetry diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index 8c3786d51a..62ac6e1ea4 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -44,7 +44,6 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG) if ((isd > oG%isc) .or. (ied < oG%ied) .or. (jsd > oG%jsc) .or. (jed > oG%jed)) & call MOM_error(FATAL, "copy_dyngrid_to_MOM_grid called with incompatible grids.") - oG%Zd_to_m = dG%Zd_to_m do i=isd,ied ; do j=jsd,jed oG%geoLonT(i,j) = dG%geoLonT(i+ido,j+jdo) oG%geoLatT(i,j) = dG%geoLatT(i+ido,j+jdo) @@ -188,7 +187,6 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG) if ((isd > dG%isc) .or. (ied < dG%ied) .or. (jsd > dG%jsc) .or. (jed > dG%jed)) & call MOM_error(FATAL, "copy_dyngrid_to_MOM_grid called with incompatible grids.") - dG%Zd_to_m = oG%Zd_to_m do i=isd,ied ; do j=jsd,jed dG%geoLonT(i,j) = oG%geoLonT(i+ido,j+jdo) dG%geoLatT(i,j) = oG%geoLatT(i+ido,j+jdo) diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 37500d31c2..51c45bc1b9 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -132,17 +132,16 @@ module MOM_dyn_horgrid real, allocatable, dimension(:,:) :: & bathyT !< Ocean bottom depth at tracer points, in depth units. - real :: Zd_to_m = 1.0 !< The conversion factor between the units of bathyT and m. logical :: bathymetry_at_vel !< If true, there are separate values for the !! basin depths at velocity points. Otherwise the effects of !! of topography are entirely determined from thickness points. real, allocatable, dimension(:,:) :: & - Dblock_u, & !< Topographic depths at u-points at which the flow is blocked, in depth units. - Dopen_u !< Topographic depths at u-points at which the flow is open at width dy_Cu, in depth units. + Dblock_u, & !< Topographic depths at u-points at which the flow is blocked, in Z. + Dopen_u !< Topographic depths at u-points at which the flow is open at width dy_Cu, in Z. real, allocatable, dimension(:,:) :: & - Dblock_v, & !< Topographic depths at v-points at which the flow is blocked, in depth units. - Dopen_v !< Topographic depths at v-points at which the flow is open at width dx_Cv, in depth units. + Dblock_v, & !< Topographic depths at v-points at which the flow is blocked, in Z. + Dopen_v !< Topographic depths at v-points at which the flow is open at width dx_Cv, in Z. real, allocatable, dimension(:,:) :: & CoriolisBu !< The Coriolis parameter at corner points, in s-1. real, allocatable, dimension(:,:) :: & @@ -160,7 +159,7 @@ module MOM_dyn_horgrid real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain real :: Rad_Earth = 6.378e6 !< The radius of the planet in meters. - real :: max_depth !< The maximum depth of the ocean in depth units (scaled by Zd_to_m). + real :: max_depth !< The maximum depth of the ocean in Z. end type dyn_horgrid_type contains @@ -287,13 +286,13 @@ subroutine rescale_dyn_horgrid_bathymetry(G, m_in_new_units) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - if (m_in_new_units == G%Zd_to_m) return + if (m_in_new_units == 1.0) return if (m_in_new_units < 0.0) & call MOM_error(FATAL, "rescale_grid_bathymetry: Negative depth units are not permitted.") if (m_in_new_units == 0.0) & call MOM_error(FATAL, "rescale_grid_bathymetry: Zero depth units are not permitted.") - rescale = G%Zd_to_m / m_in_new_units + rescale = 1.0 / m_in_new_units do j=jsd,jed ; do i=isd,ied G%bathyT(i,j) = rescale*G%bathyT(i,j) enddo ; enddo @@ -304,7 +303,6 @@ subroutine rescale_dyn_horgrid_bathymetry(G, m_in_new_units) G%Dblock_v(i,J) = rescale*G%Dblock_v(i,J) ; G%Dopen_v(i,J) = rescale*G%Dopen_v(i,J) enddo ; enddo ; endif G%max_depth = rescale*G%max_depth - G%Zd_to_m = m_in_new_units end subroutine rescale_dyn_horgrid_bathymetry From 5810cf0c7b329c6e73007c7fab77eafc706b7143 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Nov 2018 11:23:13 -0500 Subject: [PATCH 10/11] Recast internal MOM_barotropic variables into Z Recast two internal variables used in the call to set_dtbt from within barotropic_init to have units of Z instead of m and m2 Z-1 s-2 instead of m s-2, simplifying the code, and expanding dimensional consistency testing. All answers are bitwise identical in the MOM6 test cases, including rescaling Z over a large range. --- src/core/MOM_barotropic.F90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 5c6ec52fc7..7095d370b6 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2272,10 +2272,10 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) !! the effective open face areas as a !! function of barotropic flow. real, optional, intent(in) :: gtot_est !< An estimate of the total gravitational - !! acceleration, in m s-2. + !! acceleration, in m2 Z-1 s-2. real, optional, intent(in) :: SSH_add !< An additional contribution to SSH to !! provide a margin of error when - !! calculating the external wave speed, in m. + !! calculating the external wave speed, in Z. ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & @@ -2299,7 +2299,7 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) ! order 1. For stability, this may be made larger ! than physical problem would suggest. real :: add_SSH ! An additional contribution to SSH to provide a margin of error - ! when calculating the external wave speed, in m. + ! when calculating the external wave speed, in Z. real :: min_max_dt2, Idt_max2, dtbt_max logical :: use_BT_cont type(memory_size_type) :: MS @@ -2326,7 +2326,7 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) elseif (CS%Nonlinear_continuity .and. present(eta)) then call find_face_areas(Datu, Datv, G, GV, CS, MS, eta=eta, halo=0) else - call find_face_areas(Datu, Datv, G, GV, CS, MS, halo=0, add_max=add_SSH*US%m_to_Z) + call find_face_areas(Datu, Datv, G, GV, CS, MS, halo=0, add_max=add_SSH) endif det_de = 0.0 @@ -2345,8 +2345,8 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) enddo ; enddo ; enddo else do j=js,je ; do i=is,ie - gtot_E(i,j) = gtot_est * GV%H_to_m ; gtot_W(i,j) = gtot_est * GV%H_to_m - gtot_N(i,j) = gtot_est * GV%H_to_m ; gtot_S(i,j) = gtot_est * GV%H_to_m + gtot_E(i,j) = gtot_est * GV%H_to_Z ; gtot_W(i,j) = gtot_est * GV%H_to_Z + gtot_N(i,j) = gtot_est * GV%H_to_Z ; gtot_S(i,j) = gtot_est * GV%H_to_Z enddo ; enddo endif @@ -3717,9 +3717,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, character(len=40) :: mdl = "MOM_barotropic" ! This module's name. real :: Datu(SZIBS_(G),SZJ_(G)) ! Zonal open face area in H m. real :: Datv(SZI_(G),SZJBS_(G)) ! Meridional open face area in H m. - real :: gtot_estimate ! Summing GV%g_prime gives an upper-bound estimate for pbce. + real :: gtot_estimate ! Summed GV%g_prime in m2 Z-1 s-2, to give an upper-bound estimate for pbce. real :: SSH_extra ! An estimate of how much higher SSH might get, for use - ! in calculating the safe external wave speed. + ! in calculating the safe external wave speed, in Z. real :: dtbt_input, dtbt_tmp real :: wave_drag_scale ! A scaling factor for the barotropic linear wave drag ! piston velocities. @@ -3956,7 +3956,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "An estimate of how much higher SSH might get, for use \n"//& "in calculating the safe external wave speed. The \n"//& "default is the minimum of 10 m or 5% of MAXIMUM_DEPTH.", & - units="m", default=min(10.0,0.05*G%max_depth*US%Z_to_m)) + units="m", default=min(10.0,0.05*G%max_depth*US%Z_to_m), scale=US%m_to_Z) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & @@ -4141,8 +4141,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, ! Estimate the maximum stable barotropic time step. gtot_estimate = 0.0 - do k=1,G%ke ; gtot_estimate = gtot_estimate + GV%g_prime(K)*US%m_to_Z ; enddo - call set_dtbt(G, GV, US, CS, gtot_est = gtot_estimate, SSH_add = SSH_extra) + do k=1,G%ke ; gtot_estimate = gtot_estimate + GV%g_prime(K) ; enddo + call set_dtbt(G, GV, US, CS, gtot_est=gtot_estimate, SSH_add=SSH_extra) if (dtbt_input > 0.0) then CS%dtbt = dtbt_input From 7f737646b66040b5601cc6ac08480d8a598b1363 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 15 Nov 2018 18:41:13 -0500 Subject: [PATCH 11/11] +Recast ustar_gustless into Z/s Rescaled ustar_gustless and the ustar argument to get_Langmuir_number from m/s to Z/s, and added a unit_scale_type argument to set_derived_forcing fields. Also rolled some unit conversion factors into Irho0, simplifying several lines of the code in extract_IOB_stresses. All answers are bitwise identical. --- .../coupled_driver/MOM_surface_forcing.F90 | 24 +++++++++---------- .../ice_solo_driver/MOM_surface_forcing.F90 | 14 ++++++----- config_src/mct_driver/MOM_ocean_model.F90 | 2 +- .../solo_driver/MOM_surface_forcing.F90 | 2 +- src/core/MOM_forcing_type.F90 | 11 +++++---- .../vertical/MOM_CVMix_KPP.F90 | 2 +- .../vertical/MOM_energetic_PBL.F90 | 2 +- src/user/MOM_wave_interface.F90 | 18 +++++++------- 8 files changed, 39 insertions(+), 36 deletions(-) diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index 045a9dc7de..fc646a4f56 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -794,7 +794,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the !! salinity to the right time, when it is being restored. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a !! previous call to surface_forcing_init. real, dimension(SZIB_(G),SZJ_(G)), & @@ -805,7 +805,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, optional, intent(inout) :: ustar !< The surface friction velocity, in Z s-1. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(out) :: gustless_ustar !< The surface friction velocity without - !! any contributions from gustiness, in m s-1. + !! any contributions from gustiness, in Z s-1. integer, optional, intent(in) :: tau_halo !< The halo size of wind stresses to set, 0 by default. ! Local variables @@ -817,7 +817,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, real, dimension(SZIB_(G),SZJB_(G)) :: tauy_in_B ! Meridional wind stresses (in Pa) at q points real :: gustiness ! unresolved gustiness that contributes to ustar (Pa) - real :: Irho0 ! inverse of the mean density in (m^3/kg) + real :: Irho0 ! Inverse of the mean density rescaled to (Z2 m / kg) real :: taux2, tauy2 ! squared wind stresses (Pa^2) real :: tau_mag ! magnitude of the wind stress (Pa) @@ -831,7 +831,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, Isqh = G%IscB-halo ; Ieqh = G%IecB+halo ; Jsqh = G%JscB-halo ; Jeqh = G%JecB+halo i0 = is - index_bounds(1) ; j0 = js - index_bounds(3) - Irho0 = 1.0/CS%Rho0 + Irho0 = US%m_to_Z**2 / CS%Rho0 do_ustar = present(ustar) ; do_gustless = present(gustless_ustar) @@ -943,10 +943,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0)) ) & gustiness = CS%gust(i,j) endif - ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0*IOB%stress_mag(i-i0,j-j0)) + ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*IOB%stress_mag(i-i0,j-j0)) enddo ; enddo ; endif if (do_gustless) then ; do j=js,je ; do i=is,ie - gustless_ustar(i,j) = sqrt(IOB%stress_mag(i-i0,j-j0) / CS%Rho0) + gustless_ustar(i,j) = US%m_to_Z * sqrt(IOB%stress_mag(i-i0,j-j0) / CS%Rho0) !### Change to: ! gustless_ustar(i,j) = sqrt(Irho0 * IOB%stress_mag(i-i0,j-j0)) enddo ; enddo ; endif @@ -962,8 +962,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) ) if (CS%read_gust_2d) gustiness = CS%gust(i,j) endif - if (do_ustar) ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0 * tau_mag) - if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) + if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag) + if (do_gustless) gustless_ustar(i,j) = US%m_to_Z * sqrt(tau_mag / CS%Rho0) !### Change to: ! if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag) enddo ; enddo @@ -972,8 +972,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, tau_mag = G%mask2dT(i,j) * sqrt(taux_in_A(i,j)**2 + tauy_in_A(i,j)**2) gustiness = CS%gust_const if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j) - if (do_ustar) ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0 * tau_mag) - if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) + if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag) + if (do_gustless) gustless_ustar(i,j) = US%m_to_Z * sqrt(tau_mag / CS%Rho0) !### Change to: ! if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag) enddo ; enddo @@ -991,8 +991,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, gustiness = CS%gust_const if (CS%read_gust_2d) gustiness = CS%gust(i,j) - if (do_ustar) ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0 * tau_mag) - if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) + if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag) + if (do_gustless) gustless_ustar(i,j) = US%m_to_Z * sqrt(tau_mag / CS%Rho0) !### Change to: ! if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag) enddo ; enddo diff --git a/config_src/ice_solo_driver/MOM_surface_forcing.F90 b/config_src/ice_solo_driver/MOM_surface_forcing.F90 index c0b6d7dd94..301b8a9eea 100644 --- a/config_src/ice_solo_driver/MOM_surface_forcing.F90 +++ b/config_src/ice_solo_driver/MOM_surface_forcing.F90 @@ -167,15 +167,17 @@ module MOM_surface_forcing contains -subroutine set_forcing(sfc_state, forcing, fluxes, day_start, day_interval, G, CS) +subroutine set_forcing(sfc_state, forcing, fluxes, day_start, day_interval, G, US, CS) type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces - type(forcing), intent(inout) :: fluxes - type(time_type), intent(in) :: day_start - type(time_type), intent(in) :: day_interval + type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields + type(time_type), intent(in) :: day_start !< The start time of the fluxes + type(time_type), intent(in) :: day_interval !< Length of time over which these fluxes applied type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - type(surface_forcing_CS), pointer :: CS + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(surface_forcing_CS), pointer :: CS !< pointer to control struct returned by + !! a previous surface_forcing_init call ! This subroutine calls other subroutines in this file to get surface forcing fields. ! It also allocates and initializes the fields in the flux type. @@ -282,7 +284,7 @@ subroutine set_forcing(sfc_state, forcing, fluxes, day_start, day_interval, G, C ! Fields that exist in both the forcing and mech_forcing types must be copied. if (CS%variable_winds .or. CS%first_call_set_forcing) then call copy_common_forcing_fields(forces, fluxes, G) - call set_derived_forcing_fields(forces, fluxes, G, CS%Rho0) + call set_derived_forcing_fields(forces, fluxes, G, US, CS%Rho0) endif if ((CS%variable_buoyforce .or. CS%first_call_set_forcing) .and. & diff --git a/config_src/mct_driver/MOM_ocean_model.F90 b/config_src/mct_driver/MOM_ocean_model.F90 index a7b41dbe67..1a3a06b4d8 100644 --- a/config_src/mct_driver/MOM_ocean_model.F90 +++ b/config_src/mct_driver/MOM_ocean_model.F90 @@ -542,7 +542,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & #endif endif - call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%GV%Rho0) + call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%GV%Rho0) call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid) if (OS%nstep==0) then diff --git a/config_src/solo_driver/MOM_surface_forcing.F90 b/config_src/solo_driver/MOM_surface_forcing.F90 index 1a2622cfc5..d404edf9f3 100644 --- a/config_src/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/solo_driver/MOM_surface_forcing.F90 @@ -340,7 +340,7 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US ! Fields that exist in both the forcing and mech_forcing types must be copied. if (CS%variable_winds .or. CS%first_call_set_forcing) then call copy_common_forcing_fields(forces, fluxes, G) - call set_derived_forcing_fields(forces, fluxes, G, CS%Rho0) + call set_derived_forcing_fields(forces, fluxes, G, US, CS%Rho0) endif if ((CS%variable_buoyforce .or. CS%first_call_set_forcing) .and. & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 770cfc6e35..69f39c7ddb 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -47,7 +47,7 @@ module MOM_forcing_type real, pointer, dimension(:,:) :: & ustar => NULL(), & !< surface friction velocity scale (Z/s) ustar_gustless => NULL() !< surface friction velocity scale without any - !! any augmentation for gustiness (m/s) + !! any augmentation for gustiness (Z/s) ! surface buoyancy force, used when temperature is not a state variable real, pointer, dimension(:,:) :: & @@ -1961,19 +1961,20 @@ end subroutine copy_common_forcing_fields !> This subroutine calculates certain derived forcing fields based on information !! from a mech_forcing type and stores them in a (thermodynamic) forcing type. -subroutine set_derived_forcing_fields(forces, fluxes, G, Rho0) +subroutine set_derived_forcing_fields(forces, fluxes, G, US, Rho0) type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields type(ocean_grid_type), intent(in) :: G !< grid type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, intent(in) :: Rho0 !< A reference density of seawater, in kg m-3, !! as used to calculate ustar. real :: taux2, tauy2 ! Squared wind stress components, in Pa^2. - real :: Irho0 ! Inverse of the mean density in (m^3/kg) + real :: Irho0 ! Inverse of the mean density rescaled to (Z2 m / kg) integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - Irho0 = 1.0/Rho0 + Irho0 = US%m_to_Z**2 / Rho0 if (associated(forces%taux) .and. associated(forces%tauy) .and. & associated(fluxes%ustar_gustless)) then @@ -1989,7 +1990,7 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, Rho0) G%mask2dCv(i,J) * forces%tauy(i,J)**2) / & (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) - fluxes%ustar_gustless(i,j) = sqrt(sqrt(taux2 + tauy2) / Rho0) + fluxes%ustar_gustless(i,j) = US%m_to_Z * sqrt(sqrt(taux2 + tauy2) / Rho0) !### Change to: ! fluxes%ustar_gustless(i,j) = sqrt(sqrt(taux2 + tauy2) * Irho0) enddo ; enddo diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 9c5e6bb3cc..67c4173b96 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1072,7 +1072,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF endif !For now get Langmuir number based on prev. MLD (otherwise must compute 3d LA) MLD_GUESS = max( 1.*US%m_to_Z, abs(US%m_to_Z*CS%OBLdepthprev(i,j) ) ) - call get_Langmuir_Number( LA, G, GV, US, MLD_guess, surfFricVel, I, J, & + call get_Langmuir_Number( LA, G, GV, US, MLD_guess, uStar(i,j), i, j, & H=H(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES) WAVES%La_SL(i,j)=LA endif diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 0f685693c3..9e32cd2aa9 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -352,7 +352,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS real :: I_dtrho ! 1.0 / (dt * Rho0) in m3 kg-1 s-1. This is ! used convert TKE back into ustar^3. real :: U_star ! The surface friction velocity, in Z s-1. - real :: U_Star_Mean ! The surface friction without gustiness in m s-1. + real :: U_Star_Mean ! The surface friction without gustiness in Z s-1. real :: vstar ! An in-situ turbulent velocity, in m s-1. real :: Enhance_M ! An enhancement factor for vstar, based here on Langmuir impact. real :: LA ! The Langmuir number (non-dim) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 1ede95cce5..f9934615dc 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -864,7 +864,7 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, & type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, intent(in) :: i !< Meridional index of h-point integer, intent(in) :: j !< Zonal index of h-point - real, intent(in) :: ustar !< Friction velocity (m/s) + real, intent(in) :: ustar !< Friction velocity (Z/s) real, intent(in) :: HBL !< (Positive) thickness of boundary layer (Z) logical, optional, intent(in) :: Override_MA !< Override to use misalignment in LA !! calculation. This can be used if diagnostic @@ -940,7 +940,7 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, & call Get_SL_Average_Prof( GV, Dpt_LASL, H, VS_H, LA_STKy) LA_STK = sqrt(LA_STKX**2 + LA_STKY**2) elseif (WaveMethod==LF17) then - call get_StokesSL_LiFoxKemper(ustar,hbl*LA_FracHBL, GV, US, LA_STK, LA) + call get_StokesSL_LiFoxKemper(ustar, hbl*LA_FracHBL, GV, US, LA_STK, LA) endif if (.not.(WaveMethod==LF17)) then @@ -949,11 +949,11 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, & ! there is also no good reason to cap it here other then ! to prevent large enhancements in unconstrained parts of ! the curve fit parameterizations. - LA = max(WAVES%La_min,sqrt(USTAR/(LA_STK+1.e-10))) + LA = max(WAVES%La_min, sqrt(US%Z_to_m*ustar / (LA_STK+1.e-10))) endif if (Use_MA) then - WaveDirection = atan2(LA_STKy,LA_STKx) + WaveDirection = atan2(LA_STKy, LA_STKx) LA = LA / sqrt(max(1.e-8,cos( WaveDirection - ShearDirection))) endif @@ -977,7 +977,7 @@ end subroutine get_Langmuir_Number !! - BGR remove u10 input !! - BGR note: fixed parameter values should be changed to "get_params" subroutine get_StokesSL_LiFoxKemper(ustar, hbl, GV, US, UStokes_SL, LA) - real, intent(in) :: ustar !< water-side surface friction velocity (m/s) + real, intent(in) :: ustar !< water-side surface friction velocity (Z/s) real, intent(in) :: hbl !< boundary layer depth (Z) type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1001,7 +1001,7 @@ subroutine get_StokesSL_LiFoxKemper(ustar, hbl, GV, US, UStokes_SL, LA) if (ustar > 0.0) then ! Computing u10 based on u_star and COARE 3.5 relationships - call ust_2_u10_coare3p5(ustar*sqrt(GV%Rho0/1.225), u10, GV, US) + call ust_2_u10_coare3p5(US%Z_to_m*ustar*sqrt(GV%Rho0/1.225), u10, GV, US) ! surface Stokes drift UStokes = us_to_u10*u10 ! @@ -1046,13 +1046,13 @@ subroutine get_StokesSL_LiFoxKemper(ustar, hbl, GV, US, UStokes_SL, LA) sqrt( 2.0 * PI *kstar * z0) * & erfc( sqrt( 2.0 * kstar * z0 ) ) UStokes_sl = UStokes * (0.715 + r1 + r2 + r3 + r4) - LA = sqrt(ustar/UStokes_sl) + LA = sqrt(US%Z_to_m*ustar / UStokes_sl) else UStokes_sl = 0.0 LA=1.e8 endif - return -endsubroutine Get_StokesSL_LiFoxKemper + +end subroutine Get_StokesSL_LiFoxKemper !> Get SL Averaged Stokes drift from a Stokes drift Profile subroutine Get_SL_Average_Prof( GV, AvgDepth, H, Profile, Average )