diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index c8662cba15..43de125701 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -32,8 +32,6 @@ module MOM_PressureForce_Mont logical :: tides !< If true, apply tidal momentum forcing. real :: Rho0 !< The density used in the Boussinesq !! approximation [kg m-3]. - real :: Rho_atm !< The assumed atmospheric density [kg m-3]. - !! By default, Rho_atm is 0. real :: GFS_scale !< Ratio between gravity applied to top interface and the !! gravitational acceleration of the planet [nondim]. !! Usually this ratio is 1. diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 94cf169e29..09cbd14c60 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -204,7 +204,7 @@ module MOM_variables real, pointer, dimension(:,:) :: TKE_BBL => NULL() !< A term related to the bottom boundary layer source of turbulent kinetic !! energy, currently in [Z3 T-3 ~> m3 s-3], but may at some time be changed - !! to [kg Z3 m-3 T-3 ~> W m-2]. + !! to [R Z3 T-3 ~> W m-2]. real, pointer, dimension(:,:) :: & taux_shelf => NULL(), & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa]. tauy_shelf => NULL() !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa]. diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90 index 093db28c07..0608499f92 100644 --- a/src/core/MOM_verticalGrid.F90 +++ b/src/core/MOM_verticalGrid.F90 @@ -29,7 +29,7 @@ module MOM_verticalGrid real :: mks_g_Earth !< The gravitational acceleration in unscaled MKS units [m s-2]. real :: g_Earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]. real :: Rho0 !< The density used in the Boussinesq approximation or nominal - !! density used to convert depths into mass units [kg m-3]. + !! density used to convert depths into mass units [R ~> kg m-3]. ! Vertical coordinate descriptions for diagnostics and I/O character(len=40) :: zAxisUnits !< The units that vertical coordinates are written in diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 284322f072..84c4011718 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -228,6 +228,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & real :: Rcv(SZI_(G),SZJ_(G),SZK_(G)) ! Coordinate variable potential density [R ~> kg m-3]. real :: work_3d(SZI_(G),SZJ_(G),SZK_(G)) ! A 3-d temporary work array. real :: work_2d(SZI_(G),SZJ_(G)) ! A 2-d temporary work array. + real :: rho_in_situ(SZI_(G)) ! In situ density [R ~> kg m-3] ! tmp array for surface properties real :: surface_field(SZI_(G),SZJ_(G)) @@ -328,17 +329,17 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & ! diagnose thickness/volumes of grid cells [m] if (CS%id_thkcello>0 .or. CS%id_volcello>0) then if (GV%Boussinesq) then ! thkcello = h for Boussinesq - if (CS%id_thkcello > 0) then ; if (GV%H_to_m == 1.0) then + if (CS%id_thkcello > 0) then ; if (GV%H_to_Z == 1.0) then call post_data(CS%id_thkcello, h, CS%diag) else do k=1,nz; do j=js,je ; do i=is,ie - work_3d(i,j,k) = GV%H_to_m*h(i,j,k) + work_3d(i,j,k) = GV%H_to_Z*h(i,j,k) enddo ; enddo ; enddo call post_data(CS%id_thkcello, work_3d, CS%diag) endif ; endif if (CS%id_volcello > 0) then ! volcello = h*area for Boussinesq do k=1,nz; do j=js,je ; do i=is,ie - work_3d(i,j,k) = ( GV%H_to_m*h(i,j,k) ) * US%L_to_m**2*G%areaT(i,j) + work_3d(i,j,k) = ( GV%H_to_Z*h(i,j,k) ) * US%Z_to_m*US%L_to_m**2*G%areaT(i,j) enddo ; enddo ; enddo call post_data(CS%id_volcello, work_3d, CS%diag) endif @@ -357,11 +358,11 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & do i=is,ie ! Pressure for EOS at the layer center [Pa] pressure_1d(i) = pressure_1d(i) + 0.5*GV%H_to_Pa*h(i,j,k) enddo - ! Store in-situ density [kg m-3] in work_3d + ! Store in-situ density [R ~> kg m-3] in work_3d call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, & - work_3d(:,j,k), is, ie-is+1, tv%eqn_of_state) + rho_in_situ, is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) do i=is,ie ! Cell thickness = dz = dp/(g*rho) (meter); store in work_3d - work_3d(i,j,k) = (GV%H_to_kg_m2*h(i,j,k)) / work_3d(i,j,k) + work_3d(i,j,k) = (GV%H_to_RZ*h(i,j,k)) / rho_in_situ(i) enddo do i=is,ie ! Pressure for EOS at the bottom interface [Pa] pressure_1d(i) = pressure_1d(i) + 0.5*GV%H_to_Pa*h(i,j,k) @@ -371,7 +372,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_thkcello > 0) call post_data(CS%id_thkcello, work_3d, CS%diag) if (CS%id_volcello > 0) then do k=1,nz; do j=js,je ; do i=is,ie ! volcello = dp/(rho*g)*area for non-Boussinesq - work_3d(i,j,k) = US%L_to_m**2*G%areaT(i,j) * work_3d(i,j,k) + work_3d(i,j,k) = US%Z_to_m*US%L_to_m**2*G%areaT(i,j) * work_3d(i,j,k) enddo ; enddo ; enddo call post_data(CS%id_volcello, work_3d, CS%diag) endif @@ -465,7 +466,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & !$OMP parallel do default(shared) do k=1,nz ; do j=js-1,je+1 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, & - Rcv(:,j,k), is-1, ie-is+3, tv%eqn_of_state, scale=US%kg_m3_to_R) + Rcv(:,j,k), is-1, ie-is+3, tv%eqn_of_state , scale=US%kg_m3_to_R) enddo ; enddo else ! Rcv should not be used much in this case, so fill in sensible values. do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 @@ -778,7 +779,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) z_top, & ! Height of the top of a layer or the ocean [Z ~> m]. z_bot, & ! Height of the bottom of a layer (for id_mass) or the ! (positive) depth of the ocean (for id_col_ht) [Z ~> m]. - mass, & ! integrated mass of the water column [kg m-2]. For + mass, & ! integrated mass of the water column [R Z ~> kg m-2]. For ! non-Boussinesq models this is rho*dz. For Boussinesq ! models, this is either the integral of in-situ density ! (rho*dz for col_mass) or reference density (Rho_0*dz for mass_wt). @@ -788,7 +789,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) dpress, & ! Change in hydrostatic pressure across a layer [Pa]. tr_int ! vertical integral of a tracer times density, ! (Rho_0 in a Boussinesq model) [TR kg m-2]. - real :: IG_Earth ! Inverse of gravitational acceleration [s2 m-1]. + real :: IG_Earth ! Inverse of gravitational acceleration [s2 Z m-2 ~> s2 m-1]. integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -796,7 +797,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) if (CS%id_mass_wt > 0) then do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo do k=1,nz ; do j=js,je ; do i=is,ie - mass(i,j) = mass(i,j) + GV%H_to_kg_m2*h(i,j,k) + mass(i,j) = mass(i,j) + GV%H_to_RZ*h(i,j,k) enddo ; enddo ; enddo call post_data(CS%id_mass_wt, mass, CS%diag) endif @@ -804,7 +805,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) if (CS%id_temp_int > 0) then do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo do k=1,nz ; do j=js,je ; do i=is,ie - tr_int(i,j) = tr_int(i,j) + (GV%H_to_kg_m2*h(i,j,k))*tv%T(i,j,k) + tr_int(i,j) = tr_int(i,j) + (GV%H_to_RZ*h(i,j,k))*tv%T(i,j,k) enddo ; enddo ; enddo call post_data(CS%id_temp_int, tr_int, CS%diag) endif @@ -812,7 +813,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) if (CS%id_salt_int > 0) then do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo do k=1,nz ; do j=js,je ; do i=is,ie - tr_int(i,j) = tr_int(i,j) + (GV%H_to_kg_m2*h(i,j,k))*tv%S(i,j,k) + tr_int(i,j) = tr_int(i,j) + (GV%H_to_RZ*h(i,j,k))*tv%S(i,j,k) enddo ; enddo ; enddo call post_data(CS%id_salt_int, tr_int, CS%diag) endif @@ -830,7 +831,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo if (GV%Boussinesq) then if (associated(tv%eqn_of_state)) then - IG_Earth = 1.0 / GV%mks_g_Earth + IG_Earth = 1.0 / (US%Z_to_m*GV%mks_g_Earth) ! do j=js,je ; do i=is,ie ; z_bot(i,j) = -P_SURF(i,j)/GV%H_to_Pa ; enddo ; enddo do j=G%jscB,G%jecB+1 ; do i=G%iscB,G%iecB+1 z_bot(i,j) = 0.0 @@ -844,17 +845,17 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) z_top, z_bot, 0.0, US%R_to_kg_m3*GV%Rho0, GV%mks_g_Earth*US%Z_to_m, & G%HI, G%HI, tv%eqn_of_state, dpress) do j=js,je ; do i=is,ie - mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth + mass(i,j) = mass(i,j) + dpress(i,j) * US%kg_m3_to_R*IG_Earth enddo ; enddo enddo else do k=1,nz ; do j=js,je ; do i=is,ie - mass(i,j) = mass(i,j) + (GV%H_to_m*US%R_to_kg_m3*GV%Rlay(k))*h(i,j,k) + mass(i,j) = mass(i,j) + (GV%H_to_Z*GV%Rlay(k))*h(i,j,k) enddo ; enddo ; enddo endif else do k=1,nz ; do j=js,je ; do i=is,ie - mass(i,j) = mass(i,j) + GV%H_to_kg_m2*h(i,j,k) + mass(i,j) = mass(i,j) + GV%H_to_RZ*h(i,j,k) enddo ; enddo ; enddo endif if (CS%id_col_mass > 0) then @@ -866,7 +867,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) ! pbo = (mass * g) + p_surf ! where p_surf is the sea water pressure at sea water surface. do j=js,je ; do i=is,ie - btm_pres(i,j) = mass(i,j) * GV%mks_g_Earth + btm_pres(i,j) = US%RZ_to_kg_m2*mass(i,j) * GV%mks_g_Earth if (associated(p_surf)) then btm_pres(i,j) = btm_pres(i,j) + p_surf(i,j) endif @@ -1353,20 +1354,20 @@ subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dy type(tracer_registry_type), pointer :: Reg !< Pointer to the tracer registry ! Local variables - real, dimension(SZIB_(G), SZJ_(G)) :: umo2d ! Diagnostics of integrated mass transport [kg s-1] - real, dimension(SZI_(G), SZJB_(G)) :: vmo2d ! Diagnostics of integrated mass transport [kg s-1] - real, dimension(SZIB_(G), SZJ_(G), SZK_(G)) :: umo ! Diagnostics of layer mass transport [kg s-1] - real, dimension(SZI_(G), SZJB_(G), SZK_(G)) :: vmo ! Diagnostics of layer mass transport [kg s-1] + real, dimension(SZIB_(G), SZJ_(G)) :: umo2d ! Diagnostics of integrated mass transport [R Z L2 T-1 ~> kg s-1] + real, dimension(SZI_(G), SZJB_(G)) :: vmo2d ! Diagnostics of integrated mass transport [R Z L2 T-1 ~> kg s-1] + real, dimension(SZIB_(G), SZJ_(G), SZK_(G)) :: umo ! Diagnostics of layer mass transport [R Z L2 T-1 ~> kg s-1] + real, dimension(SZI_(G), SZJB_(G), SZK_(G)) :: vmo ! Diagnostics of layer mass transport [R Z L2 T-1 ~> kg s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tend ! Change in layer thickness due to dynamics ! [H s-1 ~> m s-1 or kg m-2 s-1]. real :: Idt ! The inverse of the time interval [T-1 ~> s-1] - real :: H_to_kg_m2_dt ! A conversion factor from accumulated transports to fluxes - ! [kg L-2 H-1 T-1 ~> kg m-3 s-1 or s-1]. + real :: H_to_RZ_dt ! A conversion factor from accumulated transports to fluxes + ! [R Z H-1 T-1 ~> kg m-3 s-1 or s-1]. integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Idt = 1. / dt_trans - H_to_kg_m2_dt = GV%H_to_kg_m2 * US%L_to_m**2 * US%s_to_T * Idt + H_to_RZ_dt = GV%H_to_RZ * Idt call diag_save_grids(diag) call diag_copy_storage_to_diag(diag, diag_pre_dyn) @@ -1374,28 +1375,28 @@ subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dy if (IDs%id_umo_2d > 0) then umo2d(:,:) = 0.0 do k=1,nz ; do j=js,je ; do I=is-1,ie - umo2d(I,j) = umo2d(I,j) + uhtr(I,j,k) * H_to_kg_m2_dt + umo2d(I,j) = umo2d(I,j) + uhtr(I,j,k) * H_to_RZ_dt enddo ; enddo ; enddo call post_data(IDs%id_umo_2d, umo2d, diag) endif if (IDs%id_umo > 0) then ! Convert to kg/s. do k=1,nz ; do j=js,je ; do I=is-1,ie - umo(I,j,k) = uhtr(I,j,k) * H_to_kg_m2_dt + umo(I,j,k) = uhtr(I,j,k) * H_to_RZ_dt enddo ; enddo ; enddo call post_data(IDs%id_umo, umo, diag, alt_h = diag_pre_dyn%h_state) endif if (IDs%id_vmo_2d > 0) then vmo2d(:,:) = 0.0 do k=1,nz ; do J=js-1,je ; do i=is,ie - vmo2d(i,J) = vmo2d(i,J) + vhtr(i,J,k) * H_to_kg_m2_dt + vmo2d(i,J) = vmo2d(i,J) + vhtr(i,J,k) * H_to_RZ_dt enddo ; enddo ; enddo call post_data(IDs%id_vmo_2d, vmo2d, diag) endif if (IDs%id_vmo > 0) then ! Convert to kg/s. do k=1,nz ; do J=js-1,je ; do i=is,ie - vmo(i,J,k) = vhtr(i,J,k) * H_to_kg_m2_dt + vmo(i,J,k) = vhtr(i,J,k) * H_to_RZ_dt enddo ; enddo ; enddo call post_data(IDs%id_vmo, vmo, diag, alt_h = diag_pre_dyn%h_state) endif @@ -1501,10 +1502,11 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag diag, 'Mass of liquid ocean', 'kg', standard_name='sea_water_mass') CS%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, Time, & - long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', v_extensive=.true.) + long_name = 'Cell Thickness', standard_name='cell_thickness', & + units='m', conversion=US%Z_to_m, v_extensive=.true.) CS%id_h_pre_sync = register_diag_field('ocean_model', 'h_pre_sync', diag%axesTL, Time, & - long_name = 'Cell thickness from the previous timestep', units='m', & - v_extensive=.true., conversion=GV%H_to_m) + long_name = 'Cell thickness from the previous timestep', & + units='m', conversion=GV%H_to_m, v_extensive=.true.) ! Note that CS%id_volcello would normally be registered here but because it is a "cell measure" and ! must be registered first. We earlier stored the handle of volcello but need it here for posting @@ -1707,24 +1709,24 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag endif CS%id_mass_wt = register_diag_field('ocean_model', 'mass_wt', diag%axesT1, Time, & - 'The column mass for calculating mass-weighted average properties', 'kg m-2') + 'The column mass for calculating mass-weighted average properties', 'kg m-2', conversion=US%RZ_to_kg_m2) if (use_temperature) then CS%id_temp_int = register_diag_field('ocean_model', 'temp_int', diag%axesT1, Time, & - 'Density weighted column integrated potential temperature', 'degC kg m-2', & + 'Density weighted column integrated potential temperature', 'degC kg m-2', conversion=US%RZ_to_kg_m2, & cmor_field_name='opottempmint', & cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_potential_temperature',& cmor_standard_name='Depth integrated density times potential temperature') CS%id_salt_int = register_diag_field('ocean_model', 'salt_int', diag%axesT1, Time, & - 'Density weighted column integrated salinity', 'psu kg m-2', & + 'Density weighted column integrated salinity', 'psu kg m-2', conversion=US%RZ_to_kg_m2, & cmor_field_name='somint', & cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_salinity',& cmor_standard_name='Depth integrated density times salinity') endif CS%id_col_mass = register_diag_field('ocean_model', 'col_mass', diag%axesT1, Time, & - 'The column integrated in situ density', 'kg m-2') + 'The column integrated in situ density', 'kg m-2', conversion=US%RZ_to_kg_m2) CS%id_col_ht = register_diag_field('ocean_model', 'col_height', diag%axesT1, Time, & 'The height of the water column', 'm', conversion=US%Z_to_m) @@ -1841,16 +1843,20 @@ subroutine register_transport_diags(Time, G, GV, US, IDs, diag) 'Accumulated meridional thickness fluxes to advect tracers', 'kg', & x_cell_method='sum', v_extensive=.true., conversion=H_convert*US%L_to_m**2) IDs%id_umo = register_diag_field('ocean_model', 'umo', & - diag%axesCuL, Time, 'Ocean Mass X Transport', 'kg s-1', & + diag%axesCuL, Time, 'Ocean Mass X Transport', & + 'kg s-1', conversion=US%RZ_T_to_kg_m2s*US%L_to_m**2, & standard_name='ocean_mass_x_transport', y_cell_method='sum', v_extensive=.true.) IDs%id_vmo = register_diag_field('ocean_model', 'vmo', & - diag%axesCvL, Time, 'Ocean Mass Y Transport', 'kg s-1', & + diag%axesCvL, Time, 'Ocean Mass Y Transport', & + 'kg s-1', conversion=US%RZ_T_to_kg_m2s*US%L_to_m**2, & standard_name='ocean_mass_y_transport', x_cell_method='sum', v_extensive=.true.) IDs%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', & - diag%axesCu1, Time, 'Ocean Mass X Transport Vertical Sum', 'kg s-1', & + diag%axesCu1, Time, 'Ocean Mass X Transport Vertical Sum', & + 'kg s-1', conversion=US%RZ_T_to_kg_m2s*US%L_to_m**2, & standard_name='ocean_mass_x_transport_vertical_sum', y_cell_method='sum') IDs%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d', & - diag%axesCv1, Time, 'Ocean Mass Y Transport Vertical Sum', 'kg s-1', & + diag%axesCv1, Time, 'Ocean Mass Y Transport Vertical Sum', & + 'kg s-1', conversion=US%RZ_T_to_kg_m2s*US%L_to_m**2, & standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum') IDs%id_dynamics_h = register_diag_field('ocean_model','dynamics_h', & diag%axesTl, Time, 'Change in layer thicknesses due to horizontal dynamics', & diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index c7c8d81019..775cf39c22 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -386,7 +386,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ PE_pt ! The potential energy at each point [J]. real, dimension(SZI_(G),SZJ_(G)) :: & Temp_int, Salt_int ! Layer and cell integrated heat and salt [J] and [g Salt]. - real :: HL2_to_kg ! A conversion factor form a thickness-volume to mass [kg H-1 L-2 ~> kg m-3 or 1] + real :: HL2_to_kg ! A conversion factor from a thickness-volume to mass [kg H-1 L-2 ~> kg m-3 or 1] real :: KE_scale_factor ! The combination of unit rescaling factors in the kinetic energy ! calculation [kg T2 H-1 L-2 s-2 ~> kg m-3 or nondim] real :: PE_scale_factor ! The combination of unit rescaling factors in the potential energy diff --git a/src/ice_shelf/MOM_marine_ice.F90 b/src/ice_shelf/MOM_marine_ice.F90 index 533fb5d9ec..4e3ce7401e 100644 --- a/src/ice_shelf/MOM_marine_ice.F90 +++ b/src/ice_shelf/MOM_marine_ice.F90 @@ -96,8 +96,6 @@ subroutine iceberg_forces(G, forces, use_ice_shelf, sfc_state, & forces%rigidity_ice_v(i,J) = forces%rigidity_ice_v(i,J) + kv_rho_ice * & min(forces%mass_berg(i,j), forces%mass_berg(i,j+1)) enddo ; enddo - !### This halo update may be unnecessary. Test it. -RWH - call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, G%domain, TO_ALL, CGRID_NE) end subroutine iceberg_forces diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index b2519d47ad..63461df157 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -131,8 +131,8 @@ subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file) 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 ! Local variables - real :: g_int ! Reduced gravities across the internal interfaces [m s-2]. - real :: g_fs ! Reduced gravity across the free surface [m s-2]. + real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2]. + real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. character(len=40) :: mdl = "set_coord_from_gprime" ! This subroutine's name. integer :: k, nz nz = GV%ke @@ -165,9 +165,9 @@ subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file) 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 ! Local variables - real :: g_fs ! Reduced gravity across the free surface [m s-2]. - real :: Rlay_Ref! The surface layer's target density [kg m-3]. - real :: RLay_range ! The range of densities [kg m-3]. + real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. + real :: Rlay_Ref! The surface layer's target density [R ~> kg m-3]. + real :: RLay_range ! The range of densities [R ~> kg m-3]. character(len=40) :: mdl = "set_coord_from_layer_density" ! This subroutine's name. integer :: k, nz nz = GV%ke @@ -213,8 +213,8 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state ! Local variables real :: T_ref ! Reference temperature real :: S_ref ! Reference salinity - real :: g_int ! Reduced gravities across the internal interfaces [m s-2]. - real :: g_fs ! Reduced gravity across the free surface [m s-2]. + real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2]. + real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. character(len=40) :: mdl = "set_coord_from_TS_ref" ! This subroutine's name. integer :: k, nz nz = GV%ke @@ -263,7 +263,7 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, & real, intent(in) :: P_Ref !< The coordinate-density reference pressure [Pa]. ! Local variables real, dimension(GV%ke) :: T0, S0, Pref - real :: g_fs ! Reduced gravity across the free surface [m s-2]. + real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. integer :: k, nz character(len=40) :: mdl = "set_coord_from_TS_profile" ! This subroutine's name. character(len=200) :: filename, coord_file, inputdir ! Strings for file/path @@ -318,7 +318,7 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, & ! of the range to that in the lighter part of the range. ! Setting this greater than 1 increases the resolution for ! the denser water. - real :: g_fs ! Reduced gravity across the free surface [m s-2]. + real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. real :: a1, frac_dense, k_frac integer :: k, nz, k_light character(len=40) :: mdl = "set_coord_from_TS_range" ! This subroutine's name. @@ -390,7 +390,7 @@ subroutine set_coord_from_file(Rlay, g_prime, GV, US, param_file) 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 ! Local variables - real :: g_fs ! Reduced gravity across the free surface [m s-2]. + real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. integer :: k, nz character(len=40) :: mdl = "set_coord_from_file" ! This subroutine's name. character(len=40) :: coord_var @@ -486,7 +486,7 @@ subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file) 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 ! Local variables - real :: g_fs ! Reduced gravity across the free surface [m s-2]. + real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. character(len=40) :: mdl = "set_coord_to_none" ! This subroutine's name. integer :: k, nz nz = GV%ke diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 7f420439f3..0965bc2fd8 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2009,7 +2009,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param logical :: debug = .false. ! manually set this to true for verbose output ! data arrays - real, dimension(:), allocatable :: z_edges_in, z_in + real, dimension(:), allocatable :: z_edges_in, z_in ! Interface heights [Z ~> m] real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3] real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] @@ -2326,9 +2326,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param do k=2,nz ; Rb(k) = 0.5*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo Rb(1) = 0.0 ; Rb(nz+1) = 2.0*GV%Rlay(nz) - GV%Rlay(nz-1) - zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, Rb, G%bathyT(is:ie,js:je), & - nlevs(is:ie,js:je), nkml, nkbl, min_depth, eps_z=eps_z, & - eps_rho=eps_rho) + call find_interfaces(rho_z, z_in, kd, Rb, G%bathyT, zi, G, US, & + nlevs, nkml, nkbl, min_depth, eps_z=eps_z, eps_rho=eps_rho) if (correct_thickness) then call adjustEtaToFitBathymetry(G, GV, US, zi, h) @@ -2399,7 +2398,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param if (adjust_temperature .and. .not. useALEremapping) then call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), & - US%R_to_kg_m3*GV%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, eos) + GV%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, US, eos) endif diff --git a/src/parameterizations/vertical/MOM_full_convection.F90 b/src/parameterizations/vertical/MOM_full_convection.F90 index 5fd3d67b36..daf41a1ad3 100644 --- a/src/parameterizations/vertical/MOM_full_convection.F90 +++ b/src/parameterizations/vertical/MOM_full_convection.F90 @@ -4,9 +4,10 @@ module MOM_full_convection ! This file is part of MOM6. See LICENSE.md for the license. use MOM_grid, only : ocean_grid_type +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : int_specific_vol_dp, calculate_density_derivs +use MOM_EOS, only : calculate_density_derivs implicit none ; private @@ -17,10 +18,11 @@ module MOM_full_convection contains !> Calculate new temperatures and salinities that have been subject to full convective mixing. -subroutine full_convection(G, GV, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, & +subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, & Kddt_convect, halo) 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 [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various @@ -38,8 +40,8 @@ subroutine full_convection(G, GV, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, & ! Local variables real, dimension(SZI_(G),SZK_(G)+1) :: & - drho_dT, & ! The derivatives of density with temperature and - drho_dS ! salinity [kg m-3 degC-1] and [kg m-3 ppt-1]. + dRho_dT, & ! The derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] + dRho_dS ! The derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]. real :: h_neglect, h0 ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. ! logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. @@ -107,7 +109,7 @@ subroutine full_convection(G, GV, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, & ! These would be Te_b(:,:) = tv%T(:,j,:), etc., but the values are not used Te_b(:,:) = 0.0 ; Se_b(:,:) = 0.0 - call smoothed_dRdT_dRdS(h, tv, Kddt_smooth, drho_dT, drho_dS, G, GV, j, p_surf, halo) + call smoothed_dRdT_dRdS(h, tv, Kddt_smooth, dRho_dT, dRho_dS, G, GV, US, j, p_surf, halo) do i=is,ie do_i(i) = (G%mask2dT(i,j) > 0.0) @@ -281,8 +283,8 @@ end subroutine full_convection !! above and below, including partial calculations from a tridiagonal solver. function is_unstable(dRho_dT, dRho_dS, h_a, h_b, mix_A, mix_B, T_a, T_b, S_a, S_b, & Te_aa, Te_bb, Se_aa, Se_bb, d_A, d_B) - real, intent(in) :: dRho_dT !< The derivative of in situ density with temperature [kg m-3 degC-1] - real, intent(in) :: dRho_dS !< The derivative of in situ density with salinity [kg m-3 ppt-1] + real, intent(in) :: dRho_dT !< The derivative of in situ density with temperature [R degC-1 ~> kg m-3 degC-1] + real, intent(in) :: dRho_dS !< The derivative of in situ density with salinity [R ppt-1 ~> kg m-3 ppt-1] real, intent(in) :: h_a !< The thickness of the layer above [H ~> m or kg m-2] real, intent(in) :: h_b !< The thickness of the layer below [H ~> m or kg m-2] real, intent(in) :: mix_A !< The time integrated mixing rate of the interface above [H ~> m or kg m-2] @@ -317,7 +319,7 @@ end function is_unstable !> Returns the partial derivatives of locally referenced potential density with !! temperature and salinity after the properties have been smoothed with a small !! constant diffusivity. -subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, j, p_surf, halo) +subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, halo) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & @@ -327,10 +329,11 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, j, p_surf, halo) real, intent(in) :: Kddt !< A diffusivity times a time increment [H2 ~> m2 or kg2 m-4]. real, dimension(SZI_(G),SZK_(G)+1), & intent(out) :: dR_dT !< Derivative of locally referenced - !! potential density with temperature [kg m-3 degC-1] + !! potential density with temperature [R degC-1 ~> kg m-3 degC-1] real, dimension(SZI_(G),SZK_(G)+1), & intent(out) :: dR_dS !< Derivative of locally referenced - !! potential density with salinity [kg m-3 ppt-1] + !! potential density with salinity [R degC-1 ~> kg m-3 ppt-1] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, intent(in) :: j !< The j-point to work on. real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [Pa]. integer, optional, intent(in) :: halo !< Halo width over which to compute @@ -405,7 +408,7 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, j, p_surf, halo) do i=is,ie ; pres(i) = 0.0 ; enddo endif call calculate_density_derivs(T_f(:,1), S_f(:,1), pres, dR_dT(:,1), dR_dS(:,1), & - is-G%isd+1, ie-is+1, tv%eqn_of_state) + is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*GV%H_to_Pa ; enddo do K=2,nz do i=is,ie @@ -413,11 +416,11 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, j, p_surf, halo) S_EOS(i) = 0.5*(S_f(i,k-1) + S_f(i,k)) enddo call calculate_density_derivs(T_EOS, S_EOS, pres, dR_dT(:,K), dR_dS(:,K), & - is-G%isd+1, ie-is+1, tv%eqn_of_state) + is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*GV%H_to_Pa ; enddo enddo call calculate_density_derivs(T_f(:,nz), S_f(:,nz), pres, dR_dT(:,nz+1), dR_dS(:,nz+1), & - is-G%isd+1, ie-is+1, tv%eqn_of_state) + is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) end subroutine smoothed_dRdT_dRdS diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 06b7f0f2a5..8b96c87320 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -349,7 +349,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif call cpu_clock_begin(id_clock_kappaShear) if (CS%Vertex_shear) then - call full_convection(G, GV, h, tv, T_adj, S_adj, fluxes%p_surf, & + call full_convection(G, GV, US, h, tv, T_adj, S_adj, fluxes%p_surf, & (GV%Z_to_H**2)*kappa_dt_fill, halo=1) call calc_kappa_shear_vertex(u, v, h, T_adj, S_adj, tv, fluxes%p_surf, visc%Kd_shear, & @@ -567,7 +567,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif if (CS%user_change_diff) then - call user_change_diff(h, tv, G, GV, CS%user_change_diff_CSp, Kd_lay, Kd_int, & + call user_change_diff(h, tv, G, GV, US, CS%user_change_diff_CSp, Kd_lay, Kd_int, & T_f, S_f, dd%Kd_user) endif diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 128427c683..aaa670070b 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -30,7 +30,7 @@ module MOM_tracer_Z_init 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 + 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)), & @@ -612,63 +612,64 @@ function find_limited_slope(val, e, k) result(slope) end function find_limited_slope !> Find interface positions corresponding to density profile -function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps_z, eps_rho) result(zi) - real, dimension(:,:,:), & - intent(in) :: rho !< potential density in z-space [kg m-3 or R ~> kg m-3] - real, dimension(size(rho,3)), & - intent(in) :: zin !< Input data levels [m or Z ~> m]. - real, dimension(:), intent(in) :: Rb !< target interface densities [kg m-3 or R ~> kg m-3] - real, dimension(size(rho,1),size(rho,2)), & - intent(in) :: depth !< ocean depth [Z ~> m]. - integer, dimension(size(rho,1),size(rho,2)), & - optional, intent(in) :: nlevs !< number of valid points in each column - logical, optional, intent(in) :: debug !< optional debug flag - integer, optional, intent(in) :: nkml !< number of mixed layer pieces - integer, optional, intent(in) :: nkbl !< number of buffer layer pieces - real, optional, intent(in) :: hml !< mixed layer depth [Z ~> m]. - real, optional, intent(in) :: eps_z !< A negligibly small layer thickness [m or Z ~> m]. - real, optional, intent(in) :: eps_rho !< A negligibly small density difference [kg m-3 or R ~> kg m-3]. - real, dimension(size(rho,1),size(rho,2),size(Rb,1)) :: zi !< The returned interface, in the same units az zin. +subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, US, nlevs, nkml, nkbl, hml, debug, eps_z, eps_rho) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + integer, intent(in) :: nk_data !< The number of levels in the input data + real, dimension(SZI_(G),SZJ_(G),nk_data), & + intent(in) :: rho !< Potential density in z-space [R ~> kg m-3] + real, dimension(nk_data), intent(in) :: zin !< Input data levels [Z ~> m]. + real, dimension(SZK_(G)+1), intent(in) :: Rb !< target interface densities [R ~> kg m-3] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth !< ocean depth [Z ~> m]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), & + intent(out) :: zi !< The returned interface heights [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: nlevs !< number of valid points in each column + logical, optional, intent(in) :: debug !< optional debug flag + integer, optional, intent(in) :: nkml !< number of mixed layer pieces + integer, optional, intent(in) :: nkbl !< number of buffer layer pieces + real, optional, intent(in) :: hml !< mixed layer depth [Z ~> m]. + real, optional, intent(in) :: eps_z !< A negligibly small layer thickness [Z ~> m]. + real, optional, intent(in) :: eps_rho !< A negligibly small density difference [R ~> kg m-3]. ! Local variables - real, dimension(size(rho,1),size(rho,3)) :: rho_ ! A slice of densities [R ~> kg m-3] - real, dimension(size(rho,1)) :: depth_ + real, dimension(SZI_(G),nk_data) :: rho_ ! A slice of densities [R ~> kg m-3] logical :: unstable integer :: dir - integer, dimension(size(rho,1),size(Rb,1)) :: ki_ - real, dimension(size(rho,1),size(Rb,1)) :: zi_ - integer, dimension(size(rho,1),size(rho,2)) :: nlevs_data - integer, dimension(size(rho,1)) :: lo, hi + integer, dimension(SZI_(G),SZK_(G)+1) :: ki_ + real, dimension(SZI_(G),SZK_(G)+1) :: zi_ + integer, dimension(SZI_(G),SZJ_(G)) :: nlevs_data + integer, dimension(SZI_(G)) :: lo, hi real :: slope,rsm,drhodz,hml_ - integer :: n,i,j,k,l,nx,ny,nz,nt - integer :: nlay,kk,nkml_,nkbl_ - logical :: debug_ = .false. real :: epsln_Z ! A negligibly thin layer thickness [m or Z ~> m]. real :: epsln_rho ! A negligibly small density change [kg m-3 or R ~> kg m-3]. real, parameter :: zoff=0.999 + integer :: kk,nkml_,nkbl_ + logical :: debug_ = .false. + integer :: i, j, k, m, n, is, ie, js, je, nz - nlay=size(Rb)-1 + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke zi(:,:,:) = 0.0 if (PRESENT(debug)) debug_=debug - nx = size(rho,1); ny=size(rho,2); nz = size(rho,3) - nlevs_data(:,:) = size(rho,3) + nlevs_data(:,:) = nz nkml_ = 0 ; if (PRESENT(nkml)) nkml_ = max(0, nkml) nkbl_ = 0 ; if (PRESENT(nkbl)) nkbl_ = max(0, nkbl) hml_ = 0.0 ; if (PRESENT(hml)) hml_ = hml - epsln_Z = 1.0e-10 ; if (PRESENT(eps_z)) epsln_Z = eps_z - epsln_rho = 1.0e-10 ; if (PRESENT(eps_rho)) epsln_rho = eps_rho + epsln_Z = 1.0e-10*US%m_to_Z ; if (PRESENT(eps_z)) epsln_Z = eps_z + epsln_rho = 1.0e-10*US%kg_m3_to_R ; if (PRESENT(eps_rho)) epsln_rho = eps_rho if (PRESENT(nlevs)) then nlevs_data(:,:) = nlevs(:,:) endif - do j=1,ny + do j=js,je rho_(:,:) = rho(:,j,:) - i_loop: do i=1,nx + i_loop: do i=is,ie if (debug_) then print *,'looking for interfaces, i,j,nlevs= ',i,j,nlevs_data(i,j) print *,'initial density profile= ', rho_(i,:) @@ -712,60 +713,73 @@ function find_interfaces(rho, zin, Rb, depth, nlevs, nkml, nkbl, hml, debug, eps ki_(:,:) = 0 zi_(:,:) = 0.0 - depth_(:) = -1.0*depth(:,j) lo(:) = 1 hi(:) = nlevs_data(:,j) ki_ = bisect_fast(rho_, Rb, lo, hi) ki_(:,:) = max(1, ki_(:,:)-1) - do i=1,nx - do l=2,nlay - slope = (zin(ki_(i,l)+1) - zin(ki_(i,l))) / max(rho_(i,ki_(i,l)+1) - rho_(i,ki_(i,l)),epsln_rho) - zi_(i,l) = -1.0*(zin(ki_(i,l)) + slope*(Rb(l)-rho_(i,ki_(i,l)))) - zi_(i,l) = max(zi_(i,l), depth_(i)) - zi_(i,l) = min(zi_(i,l), -1.0*hml_) + do i=is,ie + do m=2,nz + slope = (zin(ki_(i,m)+1) - zin(ki_(i,m))) / max(rho_(i,ki_(i,m)+1) - rho_(i,ki_(i,m)),epsln_rho) + zi_(i,m) = -1.0*(zin(ki_(i,m)) + slope*(Rb(m)-rho_(i,ki_(i,m)))) + zi_(i,m) = max(zi_(i,m), -depth(i,j)) + zi_(i,m) = min(zi_(i,m), -1.0*hml_) enddo - zi_(i,nlay+1) = depth_(i) - do l=2,nkml_+1 - zi_(i,l) = max(hml_*((1.0-real(l))/real(nkml_)), depth_(i)) + zi_(i,nz+1) = -depth(i,j) + do m=2,nkml_+1 + zi_(i,m) = max(hml_*((1.0-real(m))/real(nkml_)), -depth(i,j)) enddo - do l=nlay,nkml_+2,-1 - if (zi_(i,l) < zi_(i,l+1) + epsln_Z) zi_(i,l) = zi_(i,l+1) + epsln_Z - if (zi_(i,l) > -1.0*hml_) zi_(i,l) = max(-1.0*hml_, depth_(i)) + do m=nz,nkml_+2,-1 + if (zi_(i,m) < zi_(i,m+1) + epsln_Z) zi_(i,m) = zi_(i,m+1) + epsln_Z + if (zi_(i,m) > -1.0*hml_) zi_(i,m) = max(-1.0*hml_, -depth(i,j)) enddo enddo zi(:,j,:) = zi_(:,:) enddo -end function find_interfaces +end subroutine find_interfaces !> This subroutine determines the potential temperature and salinity that !! is consistent with the target density using provided initial guess -subroutine determine_temperature(temp, salt, R, p_ref, niter, land_fill, h, k_start, eos) +subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, k_start, US, eos, h_massless) real, dimension(:,:,:), intent(inout) :: temp !< potential temperature [degC] real, dimension(:,:,:), intent(inout) :: salt !< salinity [PSU] - real, dimension(size(temp,3)), intent(in) :: R !< desired potential density [kg m-3]. + real, dimension(size(temp,3)), intent(in) :: R_tgt !< desired potential density [R ~> kg m-3]. real, intent(in) :: p_ref !< reference pressure [Pa]. integer, intent(in) :: niter !< maximum number of iterations integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer) real, intent(in) :: land_fill !< land fill value - real, dimension(:,:,:), intent(in) :: h !< layer thickness, used only to avoid working on massless layers + real, dimension(:,:,:), intent(in) :: h !< layer thickness, used only to avoid working on + !! massless layers [H ~> m or kg m-2] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(eos_type), pointer :: eos !< seawater equation of state control structure + real, optional, intent(in) :: h_massless !< A threshold below which a layer is + !! determined to be massless [H ~> m or kg m-2] real, parameter :: T_max = 31.0, T_min = -2.0 ! Local variables (All of which need documentation!) - real(kind=8), dimension(size(temp,1),size(temp,3)) :: T, S, dT, dS, rho, hin - real(kind=8), dimension(size(temp,1),size(temp,3)) :: drho_dT, drho_dS - real(kind=8), dimension(size(temp,1)) :: press + real, dimension(size(temp,1),size(temp,3)) :: & + T, S, dT, dS, & + rho, & ! Layer densities [R ~> kg m-3] + hin, & ! Input layer thicknesses [H ~> m or kg m-2] + drho_dT, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] + drho_dS ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1] + real, dimension(size(temp,1)) :: press integer :: nx, ny, nz, nt, i, j, k, n, itt real :: dT_dS_gauge ! The relative penalizing of temperature to salinity changes when ! minimizing property changes while correcting density [degC ppt-1]. real :: I_denom ! The inverse of the magnitude squared of the density gradient in - ! T-S space streched with dT_dS_gauge [m6 kg-2 ppt-1] + ! T-S space streched with dT_dS_gauge [ppt2 R-2 ~> ppt2 m6 kg-2] logical :: adjust_salt, old_fit - real, parameter :: S_min = 0.5, S_max=65.0 - real, parameter :: tol_T=1.e-4, tol_S=1.e-4, tol_rho=1.e-4 - real, parameter :: max_t_adj=1.0, max_s_adj = 0.5 - + real :: S_min, S_max + real :: tol_T ! The tolerance for temperature matches [degC] + real :: tol_S ! The tolerance for salinity matches [ppt] + real :: tol_rho ! The tolerance for density matches [R ~> kg m-3] + real :: max_t_adj, max_s_adj + + ! These hard coded parameters need to be set properly. + S_min = 0.5 ; S_max = 65.0 + max_t_adj = 1.0 ; max_s_adj = 0.5 + tol_T=1.e-4 ; tol_S=1.e-4 ; tol_rho = 1.e-4*US%kg_m3_to_R old_fit = .true. ! reproduces siena behavior ! ### The whole determine_temperature subroutine needs to be reexamined, both the algorithms @@ -780,28 +794,29 @@ subroutine determine_temperature(temp, salt, R, p_ref, niter, land_fill, h, k_st do j=1,ny dS(:,:) = 0. ! Needs to be zero everywhere since there is a maxval(abs(dS)) later... - T=temp(:,j,:) - S=salt(:,j,:) - hin=h(:,j,:) - dT=0.0 + T(:,:) = temp(:,j,:) + S(:,:) = salt(:,j,:) + hin(:,:) = h(:,j,:) + dT(:,:) = 0.0 adjust_salt = .true. iter_loop: do itt = 1,niter do k=1, nz - call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos) - call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, eos) + call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, scale=US%kg_m3_to_R) + call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, & + eos, scale=US%kg_m3_to_R) enddo do k=k_start,nz ; do i=1,nx -! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln) then - if (abs(rho(i,k)-R(k))>tol_rho) then +! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln) then + if (abs(rho(i,k)-R_tgt(k))>tol_rho) then if (old_fit) then - dT(i,k) = max(min((R(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj) + dT(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj) T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) else dT_dS_gauge = 10.0 ! 10 degC is weighted equivalently to 1 ppt. I_denom = 1.0 / (drho_dS(i,k)**2 + dT_dS_gauge**2*drho_dT(i,k)**2) - dS(i,k) = (R(k)-rho(i,k)) * drho_dS(i,k) * I_denom - dT(i,k) = (R(k)-rho(i,k)) * dT_dS_gauge**2*drho_dT(i,k) * I_denom + dS(i,k) = (R_tgt(k)-rho(i,k)) * drho_dS(i,k) * I_denom + dT(i,k) = (R_tgt(k)-rho(i,k)) * dT_dS_gauge**2*drho_dT(i,k) * I_denom T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) @@ -816,21 +831,22 @@ subroutine determine_temperature(temp, salt, R, p_ref, niter, land_fill, h, k_st if (adjust_salt .and. old_fit) then ; do itt = 1,niter do k=1, nz - call calculate_density(T(:,k),S(:,k),press,rho(:,k),1,nx,eos) - call calculate_density_derivs(T(:,k),S(:,k),press,drho_dT(:,k),drho_dS(:,k),1,nx,eos) + call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, scale=US%kg_m3_to_R) + call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, & + eos, scale=US%kg_m3_to_R) enddo do k=k_start,nz ; do i=1,nx -! if (abs(rho(i,k)-R(k))>tol .and. hin(i,k)>epsln .and. abs(T(i,k)-land_fill) < epsln ) then - if (abs(rho(i,k)-R(k)) > tol_rho) then - dS(i,k) = max(min((R(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj) +! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then + if (abs(rho(i,k)-R_tgt(k)) > tol_rho) then + dS(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj) S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) endif enddo ; enddo if (maxval(abs(dS)) < tol_S) exit enddo ; endif - temp(:,j,:)=T(:,:) - salt(:,j,:)=S(:,:) + temp(:,j,:) = T(:,:) + salt(:,j,:) = S(:,:) enddo end subroutine determine_temperature diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index ec7f907fd1..70d89497da 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -98,7 +98,7 @@ subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) ! Set whichever fluxes are to be used here. Any fluxes that ! are always zero do not need to be changed here. do j=js,je ; do i=is,ie - ! Fluxes of fresh water through the surface are in units of [kg m-2 s-1] + ! Fluxes of fresh water through the surface are in units of [R Z T-1 ~> kg m-2 s-1] ! and are positive downward - i.e. evaporation should be negative. fluxes%evap(i,j) = -0.0 * G%mask2dT(i,j) fluxes%lprec(i,j) = 0.0 * G%mask2dT(i,j) @@ -151,7 +151,7 @@ subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) Temp_restore = 0.0 do j=js,je ; do i=is,ie ! Set density_restore to an expression for the surface potential - ! density [kg m-3] that is being restored toward. + ! density [R ~> kg m-3] that is being restored toward. if (G%geoLatT(i,j) < CS%lfrslat) then Temp_restore = CS%SST_s elseif (G%geoLatT(i,j) > CS%lfrnlat) then diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 642ed41d88..6d307f843a 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -365,7 +365,6 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC ! Local variables real :: T(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for temp [degC] real :: S(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for salt [ppt] - real :: RHO(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for RHO [kg m-3] real :: h(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for thickness [H ~> m or kg m-2]. real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! A temporary array for thickness [Z ~> m] real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1]. diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index 5cd75725e3..dd7309265f 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -351,13 +351,11 @@ end subroutine Phillips_initialize_topography !! The one argument passed to initialize, Time, is set to the !! current time of the simulation. The fields which are initialized !! here are: -!! u - Zonal velocity [m s-1]. -!! v - Meridional velocity [m s-1]. -!! h - Layer thickness in m. (Must be positive.) -!! D - Basin depth in m. (Must be positive.) -!! f - The Coriolis parameter [s-1]. -!! g - The reduced gravity at each interface [m s-2] -!! Rlay - Layer potential density (coordinate variable) [kg m-3]. +!! u - Zonal velocity [L T-1 ~> m s-1]. +!! v - Meridional velocity [L T-1 ~> m s-1]. +!! h - Layer thickness [H ~> m or kg m-2] (must be positive) +!! D - Basin depth [Z ~> m] (positive downward) +!! f - The Coriolis parameter [T-1 ~> s-1]. !! If ENABLE_THERMODYNAMICS is defined: !! T - Temperature [degC]. !! S - Salinity [ppt]. diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index 3478415c60..5641035ded 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -232,7 +232,7 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, & !! only read parameters without changing h. ! Local variables real :: T0(SZK_(G)), S0(SZK_(G)) - real :: pres(SZK_(G)) ! Reference pressure [kg m-3]. + real :: pres(SZK_(G)) ! Reference pressure [Pa]. real :: drho_dT(SZK_(G)) ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]. real :: drho_dS(SZK_(G)) ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]. real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 [R ~> kg m-3]. diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index 63f8009235..c1f615fe2a 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -63,9 +63,6 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) ! Local variables real :: Temp_restore ! The temperature that is being restored toward [degC]. real :: Salin_restore ! The salinity that is being restored toward [ppt]. - real :: density_restore ! The potential density that is being restored - ! toward [kg m-3]. - real :: rhoXcp ! The mean density times the heat capacity [J m-3 degC-1]. integer :: i, j, is, ie, js, je integer :: isd, ied, jsd, jed @@ -97,7 +94,7 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) ! Set whichever fluxes are to be used here. Any fluxes that ! are always zero do not need to be changed here. do j=js,je ; do i=is,ie - ! Fluxes of fresh water through the surface are in units of [kg m-2 s-1] + ! Fluxes of fresh water through the surface are in units of [R Z T-1 ~> kg m-2 s-1] ! and are positive downward - i.e. evaporation should be negative. fluxes%evap(i,j) = -0.0 * G%mask2dT(i,j) fluxes%lprec(i,j) = 0.0 * G%mask2dT(i,j) @@ -121,8 +118,6 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) if (CS%use_temperature .and. CS%restorebuoy) then do j=js,je ; do i=is,ie - ! Set density_restore to an expression for the surface potential - ! density [kg m-3] that is being restored toward. if (CS%forcing_mask(i,j)>0.) then fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)) * & ((CS%S_restore(i,j) - state%SSS(i,j)) / (0.5 * (CS%S_restore(i,j) + state%SSS(i,j)))) diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90 index 10d04af0c3..86f3e6e99a 100644 --- a/src/user/user_change_diffusivity.F90 +++ b/src/user/user_change_diffusivity.F90 @@ -31,7 +31,7 @@ module user_change_diffusivity !! a diffusivity scaled by Kd_add is added [degLat]. real :: rho_range(4) !< 4 values that define the coordinate potential !! density range over which a diffusivity scaled by - !! Kd_add is added [kg m-3]. + !! Kd_add is added [R ~> kg m-3]. logical :: use_abs_lat !< If true, use the absolute value of latitude when !! setting lat_range. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to @@ -44,13 +44,14 @@ module user_change_diffusivity !! main code to alter the diffusivities as needed. The specific example !! implemented here augments the diffusivity for a specified range of latitude !! and coordinate potential density. -subroutine user_change_diff(h, tv, G, GV, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_add) +subroutine user_change_diff(h, tv, G, GV, US, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_add) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers !! to any available thermodynamic !! fields. Absent fields have NULL ptrs. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(user_change_diff_CS), pointer :: CS !< This module's control structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity of !! each layer [Z2 T-1 ~> m2 s-1]. @@ -64,7 +65,7 @@ subroutine user_change_diff(h, tv, G, GV, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_a !! diffusivity that is being added at !! each interface [Z2 T-1 ~> m2 s-1]. ! Local variables - real :: Rcv(SZI_(G),SZK_(G)) ! The coordinate density in layers [kg m-3]. + real :: Rcv(SZI_(G),SZK_(G)) ! The coordinate density in layers [R ~> kg m-3]. real :: p_ref(SZI_(G)) ! An array of tv%P_Ref pressures. real :: rho_fn ! The density dependence of the input function, 0-1 [nondim]. real :: lat_fn ! The latitude dependence of the input function, 0-1 [nondim]. @@ -106,13 +107,13 @@ subroutine user_change_diff(h, tv, G, GV, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_a do j=js,je if (present(T_f) .and. present(S_f)) then do k=1,nz - call calculate_density(T_f(:,j,k),S_f(:,j,k),p_ref,Rcv(:,k),& - is,ie-is+1,tv%eqn_of_state) + call calculate_density(T_f(:,j,k), S_f(:,j,k), p_ref, Rcv(:,k),& + is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) enddo else do k=1,nz - call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p_ref,Rcv(:,k),& - is,ie-is+1,tv%eqn_of_state) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rcv(:,k),& + is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) enddo endif @@ -135,7 +136,6 @@ subroutine user_change_diff(h, tv, G, GV, CS, Kd_lay, Kd_int, T_f, S_f, Kd_int_a else lat_fn = val_weights(G%geoLatT(i,j), CS%lat_range) endif - ! rho_int = 0.5*(Rcv(i,k-1) + Rcv(i,k)) rho_fn = val_weights( 0.5*(Rcv(i,k-1) + Rcv(i,k)), CS%rho_range) if (rho_fn * lat_fn > 0.0) then Kd_int(i,j,K) = Kd_int(i,j,K) + CS%Kd_add * rho_fn * lat_fn @@ -163,9 +163,9 @@ end function range_OK !! hit 0 and 1. The values in range must be in ascending order, as can be !! checked by calling range_OK. function val_weights(val, range) result(ans) - real, intent(in) :: val !< Value for which we need an answer. - real, dimension(4), intent(in) :: range !< Range over which the answer is non-zero. - real :: ans !< Return value. + real, intent(in) :: val !< Value for which we need an answer [arbitrary units]. + real, dimension(4), intent(in) :: range !< Range over which the answer is non-zero [arbitrary units]. + real :: ans !< Return value [nondim]. ! Local variables real :: x ! A nondimensional number between 0 and 1. @@ -238,7 +238,7 @@ subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS) "is applied. The four values specify the density at "//& "which the extra diffusivity starts to increase from 0, "//& "hits its full value, starts to decrease again, and is "//& - "back to 0.", units="kg m-3", default=-1.0e9) + "back to 0.", units="kg m-3", default=-1.0e9, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "USER_KD_ADD_USE_ABS_LAT", CS%use_abs_lat, & "If true, use the absolute value of latitude when "//& "checking whether a point fits into range of latitudes.", & diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index 7db78f2454..55c609802e 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -242,10 +242,10 @@ end subroutine write_user_log !! !! This subroutine initializes the fields for the simulations. !! The one argument passed to initialize, Time, is set to the -!! current time of the simulation. The fields which are initialized +!! current time of the simulation. The fields which might be initialized !! here are: -!! - u - Zonal velocity [m s-1]. -!! - v - Meridional velocity [m s-1]. +!! - u - Zonal velocity [Z T-1 ~> m s-1]. +!! - v - Meridional velocity [Z T-1 ~> m s-1]. !! - h - Layer thickness [H ~> m or kg m-2]. (Must be positive.) !! - G%bathyT - Basin depth [Z ~> m]. (Must be positive.) !! - G%CoriolisBu - The Coriolis parameter [T-1 ~> s-1]. @@ -255,7 +255,7 @@ end subroutine write_user_log !! - T - Temperature [degC]. !! - S - Salinity [psu]. !! If BULKMIXEDLAYER is defined: -!! - Rml - Mixed layer and buffer layer potential densities [kg m-3]. +!! - Rml - Mixed layer and buffer layer potential densities [R ~> kg m-3]. !! If SPONGE is defined: !! - A series of subroutine calls are made to set up the damping !! rates and reference profiles for all variables that are damped