module ocean_density_mod #define COMP isc:iec,jsc:jec ! ! S.M. Griffies ! ! ! Russell Fiedler ! ! ! ! Compute density and related quantities. ! ! ! ! ! This module computes the in-situ density and its partial derivatives with ! respect to conservative temperature or potential temperature, and with ! respect to salinity. ! ! There are three basic means for performing this calculation. ! ! A/ Linear equation for use in idealized studies ! ! This equation renders density a linear function of potential ! temperature and salinity. All nonlinearities are ignored, as are ! pressure effects. ! ! The valid range for theta and salinity arbitrary for the ! linear equation of state. ! ! B/ pre-TEOS10 method: this method uses density as a rational ! polynomial function of potential temperature, practical salinity, ! and gauge pressure. There is also an implementation that computes ! density as a function of conservative temperature rather than ! potential temperature. The equation of state is based on that ! documented in Jackett, McDougall, Feistel, Wright, and Griffies(2006). ! ! This equation of state is valid over a "cone-shaped" range ! corresponding to ! ! 0psu <= salinity <= 40 psu ! ! -3C <= theta <= 40C "theta" = either conservative or potential temp ! ! 0dbar <= pressure <= 8000dbar ! ! with the cone getting smaller in the deeper ocean where ! theta and salinity vary over a smaller range. ! ! Input variables are the following: ! ! --salinity in psu or g/kg ! --conservative temperature or potential temperature (theta) in deg C ! --pressure in dbars (1bar = 10dbar = 10^5 Newton/m^2 = 10^5 Pascals). ! ! Note that in the ocean, pressure increases roughly by 1dbar for each meter depth. ! Also note that pressure is the "sea pressure", which is the absolute pressure ! minus the pressure of a standard atmosphere, which is 10.1325 dbars. ! ! check values
! ! for "theta" = conservative temperature ! rho(s=20psu,theta=20C,p=1000dbar) = 1017.842890411975 (kg/m^3)
! alpha(s=20psu,theta=20C,p=1000dbar) = 2.436057013634649e-4 (1/C)
! beta(s=20psu,theta=20C,p=1000dbar) = 7.314818108935248e-4 (1/psu)
! ! for "theta" = potential temperature ! rho(s=20psu,theta=20C,p=1000dbar) = 1017.728868019642 (kg/m^3)
! alpha(s=20psu,theta=20C,p=1000dbar) = 2.525481286927133e-4 (1/C)
! beta(s=20psu,theta=20C,p=1000dbar) = 7.379638527217575e-4 (1/psu)
! ! This equation of state should be suitable for purposes of realistic ! global ocean climate modeling. ! ! C/ TEOS10 method: this method makes use of the recommendations from ! the SCOR working group on seawater thermodynamics, 2010. Here, density ! is a function of conservative temperature and absolute salinity. ! The equation is valid from 0 g/kg salinty to a very large value. ! ! ! ! ! ! ! ! McDougall T.J., P.M. Barker, R. Feistel and D.R. Jackett, 2011: A ! computationally efficient 48-term expression for the density of ! seawater in terms of Conservative Temperature, and related properties ! of seawater. To be submitted to Ocean Science. ! ! ! ! IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of ! seawater 2010: Calculation and use of thermodynamic properties. ! Intergovernmental Oceanographic Commission, Manuals and Guides No. ! 56, UNESCO (English), 196 pp. ! ! ! ! D. Iudicone, G. Madec, and T.J. McDougall (2008) ! Water-mass transformations in a neutral density framework and the ! key role of light penetration. JPO vol 38, pages 1357-1376. ! ! ! ! Jackett, McDougall, Feistel, Wright, and Griffies (2006) ! Algorithms for density, potential temperature, conservative ! temperature, and freezing temperature of seawater. ! Journal of Atmospheric and Oceanic Technology, 2006, ! volume 23, pages 1709-1728. ! ! ! ! McDougall and Jackett (2005) ! The material derivative of neutral density ! Journal of Marine Research, vol 63, pages 159-185. ! ! ! ! Feistel (2003), A new extended Gibbs thermodynamic potential ! of seawater. Progress in Oceanography. vol 58, pages 43-114. ! ! ! ! S.M. Griffies, M.J. Harrison, R.C. Pacanowski, and A. Rosati ! A Technical Guide to MOM4 (2003) ! ! ! ! S.M. Griffies, R.C. Pacanowski, R.M. Schmidt, and V. Balaji ! Tracer Conservation with an Explicit Free Surface Method for ! Z-coordinate Ocean Models ! Monthly Weather Review (2001) vol 129 pages 1081--1098 ! ! ! ! T. McDougall (1987) ! Cabbeling, Thermobaricity, and water mass conversion ! JGR vol 92, pages 5448-5464 ! ! ! ! ! Density is computed as a function of conservative temperature (degC) ! or potential temperature (degC), salinity (psu or g/kg), and pressure (dbar). ! The pressure contribution includes that from the free surface height ! and the applied atmospheric and/or sea ice pressure. However, it is referenced ! to standard atmosphere, so that we use the "gauge" pressure rather than the ! full in-situ pressure. ! ! For vert_coordinate==GEOPOTENTIAL, ZSTAR, or ZSIGMA, baroclinic component of ! hydrostatic pressure is not known until density is known. In this case, ! the baroclinic pressure contribution to density is lagged by a time step. ! rho(tau) = rho[theta(tau),s(tau), p_atm(tau) + p_fs(tau) + p_baroclinic(tau-1)]. ! This issue does not arise when using vert_coordinate=PRESSURE, PSTAR, or PSIGMA. ! ! ! ! ! ! ! ! ! Set true to write a restart. False setting only for rare ! cases where wish to benchmark model without measuring the cost ! of writing restarts and associated chksums. ! Default is write_a_restart=.true. ! ! ! ! Standard atmospheric pressure (dbar). The realistic ! EOS used in MOM requires "sea pressure" as an argument ! rather than absolute pressure. Sea pressure is ! absolute pressure minus a standard atmospheric pressure ! of 10.1325dbar. ! ! For models that do have a realistic atmospheric loading, then it ! is appropriate to remove 10.1325dbar prior to computing the EOS. ! For those cases with zero atmospheric pressure, then it is not ! necessary to remove the standard atmosphere. The default for the ! press_standard is 0.0dbar. ! ! ! ! Conservative temperature or potential temperature for ! testing the EOS. ! ! ! Salinity for testing the EOS. ! ! ! Sea pressure for testing the EOS. ! ! ! ! Conservative temperature or potential temperature for ! testing the equation for neutral density. ! ! ! Salinity the equation for neutral density. ! ! ! ! Set to true to use TEOS-10 equation of state, which ! is a function of conservative temperature and absolute ! salinity. ! Default eos_teos10=.false. ! ! ! Set to true to use pre-TEOS-10 equation of state, which ! is a function of potential temperature and practical salinity, ! or conservative temperature and practical salinity. ! Default eos_preteos10=.false. ! ! ! Set to true to use an idealized linear equation of state, which ! has no pressure dependence, and is a linear function of salinity ! and temperature. ! Default eos_linear=.false. ! ! ! Constant "thermal expansion coefficient" for linear EOS ! rho = rho0 - alpha_linear_eos*theta + beta_linear_eos*salinity ! ! ! Constant "saline contraction coefficient" for linear EOS ! rho = rho0 - alpha_linear_eos*theta + beta_linear_eos*salinity ! ! ! ! Reference sea pressure for computing diagnostic potential density ! of use for computing diagnostics with potential density. ! Default potrho_press=2000.0 ! ! ! Minimum potential density used to partition vertical according ! to potential density. ! ! ! Maximum potential density used to partition vertical according ! to potential density. ! ! ! ! Set to true to compute the neutral density according to ! the omega method based on Klocker and McDougall. ! This approach has not yet been coded. Presently as a ! placeholder we use potential density referenced to 2000dbar. ! Default neutral_density_omega=.false. ! ! ! Set to true to compute the neutral density as just ! a selected potential density, set according to potrho_press. ! Since the neutral_density_omega approach has yet to be coded, ! we only have the neutral_density_potrho option to choose from ! at this time. ! Default neutral_density_potrho=.true. ! ! ! If neutral_density_potrho=true., use a different reference sea pressure ! nrho_press, instead of potrho_press? Allows multiple diagnostic potential ! densities to be output. ! Default diff_nrho_press=.false. ! ! Reference sea pressure for computing diagnostic potential density ! for the 'neutral' density variable used in computing diagnostics with nrho. ! Default nrho_press=2000.0 ! ! ! Set to true to use temperature instead of neutral density as the ! binning variable for water-mass diagnostics. ! Default neutral_density_theta=.false. ! ! ! Minimum neutral density used to partition vertical according ! to rational polynomial approximation to neutral density. ! ! ! Maximum neutral density used to partition vertical according ! to rational polynomial approximation to neutral density. ! ! ! ! Minimum conservative temperature or potential temperature used to ! partition vertical according to temperature. ! ! ! Maximum conservative temperature or potential temperature used to ! partition vertical according to temperature. ! ! ! ! Number of classes used to partition vertical according to potential density, ! conservative temperature, or potential temperature. Used for diagnostics. ! ! ! ! To smooth the vertical temp and salt derivative for diagnosing ! the buoyancy frequency. Default buoyfreq_smooth_vert=.true. ! ! ! ! To normalize the inverse vertical derivative of ! locally referenced potential density for computing the ! buoyancy frequency. Default epsln_drhodz=1e-10. ! ! ! ! To normalize the inverse vertical derivative of locally referenced ! potential density for computing neutral_rho and wdian diagnostics. ! Default epsln_drhodz_diag=1e-10. ! ! ! ! A diagnostic maximum neutral slope for use in computing which direction ! is deemed the most stratified. For use in computing the stratification_factor ! which is then used to diagnose the dianeutral mass transport. ! smax_diag should corresond to the choice used in neutral diffusion scheme. ! Should have 0 <= smax_diag <= 1.0. ! Default smax_diag=-1.0, in which case we compute the smax according to ! the vertical to horizontal grid aspect ratio. This method ensures that ! the slope is adequately "resolved" by the grid. ! ! ! To compute the diagnostic maximum neutral slope within a column as the minimum ! vertical to horizontal grid aspect ratio. This method ensures that ! the slope is adequately "resolved" by the grid, and that all depths use the ! same definition of "resolved", even if presumably thicker grid cells can ! "resolve" larger neutral slopes. This approach is not very useful generally, ! so it is retained only for testing purposes. ! Default smax_min_in_column=.false. ! ! ! ! For cases where use the domain masking, it is necessary to initialize the field ! denominator_r to nonzero in order to avoid NaNs in the case when change processor ! layout in between restarts. Note that when use solid wall boundary conditions, ! this logical should remain false in order to bitwise reproduce across restarts. ! Default mask_domain_restart=.false. ! ! ! ! When computing drhodz_diag, we can enforce that it is negative, ! thus reflecting a stable stratification. The field drhodz_diag ! is used for many water mass transformation diagnostics, such as ! wdian_rho. Allowing for unstable profiles can bias the wdian_rho ! calculation in an improper way, since the magnitude of drhodz_diag ! is very small when it is positive, whereas it is larger magnitude when ! negative. Default drhodz_diag_stable=.true. ! ! ! ! To perform the diagnostic calculation of grad_nrho_lrpotrho ! for analysis diagnostics. This factor is not well constrained, ! and can be problematic in certain regions. So presently we do ! not recommend computing it, so that the default ! is grad_nrho_lrpotrho_compute=.false. ! ! ! Maximum value used for grad_nrho_lrpotrho. ! Default grad_nrho_lrpotrho_max=10. ! ! ! Minimum value used for grad_nrho_lrpotrho. ! Default grad_nrho_lrpotrho_min=1. ! ! ! ! For doing an S2D smoothing of the stratification factor used ! for diagnostic purposes. Requires an extra call to mpp update. ! Default smooth_stratification_factor=.false. since the smoothing ! incurs a cost that should be borne only when desired. ! ! ! ! To update the watermass_factor and stratification_factor ! for use in the water mass transformation diagnostics. ! Default update_diagnostic_factors=.false. ! ! ! ! For debugging nonlinear equation of state ! ! ! ! For debugging, it is often useful to have rho=rho0 uniform. ! ! ! ! For idealized tests, set the in situ density equal to the ! potential density referenced to potrho_press. All density ! derivatives will also be computed with respect to constant ! potrho_press pressure. ! Default density_equal_potrho=.false. ! ! ! ! Set true to do bitwise exact global sum. When it is false, the global ! sum will be non-bitwise_exact, but will significantly increase ! efficiency. ! default: do_bitwise_exact_sum=.false. ! ! ! #include use constants_mod, only: epsln, c2dbars use diag_manager_mod, only: register_diag_field, register_static_field, diag_axis_init use diag_manager_mod, only: need_data, send_data use fms_mod, only: write_version_number, mpp_error use fms_mod, only: field_exist, FATAL, WARNING use fms_mod, only: open_namelist_file, check_nml_error, close_file, file_exist use fms_io_mod, only: register_restart_field, save_restart, restore_state use fms_io_mod, only: reset_field_pointer, restart_file_type use mpp_domains_mod, only: mpp_update_domains, mpp_global_sum use mpp_domains_mod, only: BITWISE_EXACT_SUM, NON_BITWISE_EXACT_SUM use mpp_mod, only: input_nml_file, stdout, stdlog use platform_mod, only: i8_kind use time_manager_mod, only: time_type, increment_time use field_manager_mod, only: fm_get_index use ocean_domains_mod, only: get_local_indices use ocean_operators_mod, only: FDX_T, FDY_T, FMX, FMY, S2D use ocean_parameters_mod, only: GEOPOTENTIAL, ZSTAR, ZSIGMA, PRESSURE, PSTAR, PSIGMA use ocean_parameters_mod, only: CONSERVATIVE_TEMP, POTENTIAL_TEMP, PRACTICAL_SALT, PREFORMED_SALT use ocean_parameters_mod, only: missing_value, onefourth, onehalf, rho0r, rho0, grav use ocean_pressure_mod, only: pressure_in_dbars, pressure_in_dbars_blob use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_thickness_type use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type, ocean_options_type use ocean_types_mod, only: ocean_prog_tracer_type, ocean_density_type, ocean_external_mode_type use ocean_types_mod, only: ocean_diag_tracer_type, ocean_lagrangian_type use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, write_chksum_3d use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk6 use ocean_workspace_mod, only: wrk1_v, wrk2_v, wrk3_v, wrk4_v use ocean_workspace_mod, only: wrk1_2d, wrk2_2d, wrk3_2d, wrk4_2d implicit none private real, dimension(:), allocatable :: a, b ! polynomial coefficients in the realistic EOS real, dimension(:,:), allocatable :: rhodz_tau ! vertical integral of rho*dzt at time tau real, dimension(:,:), allocatable :: rhodz_taup1 ! vertical integral of rho*dzt at time taup1 real, dimension(:,:,:), allocatable :: denominator_r ! reciprocal of denominator for relastic EOS integer :: vert_coordinate ! for setting which vertical coordinate to use integer :: global_sum_flag ! for determing type of global sum BITWISE/NON ! polynomial coefficients in the preTEOS10 EOS real :: a0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11 real :: b0, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11, b12 ! multiplied constants for density partial derivs in the preTEOS10 EOS real :: two_a2, three_a3, two_a6, six_a3, two_a8, two_a10, two_a11, four_a11 real :: two_b2, three_b3, six_b3, four_b4, twelve_b4, three_b7, six_b7 real :: onep5_b8, onep5_b9, two_b9, three_b9, two_b11, three_b11, six_b11, three_b12 ! polynomial coefficients in the TEOS10 eos real :: v01,v02,v03,v04,v05,v06,v07,v08,v09,v10 real :: v11,v12,v13,v14,v15,v16,v17,v18,v19,v20 real :: v21,v22,v23,v24,v25,v26,v27,v28,v29,v30 real :: v31,v32,v33,v34,v35,v36,v37,v38,v39,v40 real :: v41,v42,v43,v44,v45,v46,v47,v48 ! multiplied constants for density partial derivs in the TEOS10 EOS real :: two_v03, two_v07, two_v10, two_v14, two_v17, two_v19, two_v18, two_v20, two_v23, two_v28, two_v33 real :: two_v36, two_v39, two_v43, two_v44, two_v45, two_v46 real :: three_v04, three_v11, three_v24, three_v29, three_v34, three_v40, three_v47, three_v48 real :: four_v25, four_v30, four_v35 real :: onep5_v08, onep5_v09, onep5_v10, onep5_v11, onep5_v31, onep5_v32, onep5_v33, onep5_v34, onep5_v35 real :: six_v04, six_v11, six_v24, six_v29, six_v34, six_v40, six_v47, six_v48 real :: twelve_v25, twelve_v30, twelve_v35 real :: p75_v08, p75_v09, p75_v10, p75_v11, p75_v31, p75_v32, p75_v33, p75_v34, p75_v35 real :: three_v10, three_v33 real :: fourp5_v11, fourp5_v34 real :: four_v19, four_v45 real :: six_v35 ! polynomial coefficients in neutral density real :: a0n, a1n, a2n, a3n, a4n, a5n, a6n real :: b0n, b1n, b2n, b3n, b4n, b5n, b6n, b7n, b8n, b9n ! time steps real :: dtts=0.0 ! inverse area (1/m) of k=1 tracer cells real :: area_total_r ! inverse neutralrho_interval real :: neutralrho_interval_r ! inverse gravity constant real :: grav_r ! for diagnostics manager logical :: used integer :: id_eos_salinity =-1 integer :: id_drhodtheta =-1 integer :: id_drhodsalt =-1 integer :: id_dpotrhodtheta =-1 integer :: id_dpotrhodsalt =-1 integer :: id_dpotrhodpress =-1 integer :: id_drhodpress =-1 integer :: id_thermal_expand =-1 integer :: id_haline_contract=-1 integer :: id_sound_speed2 =-1 integer :: id_press =-1 integer :: id_rho =-1 integer :: id_rhoE =-1 integer :: id_rho0 =-1 integer :: id_neutral_rho =-1 integer :: id_pot_rho =-1 integer :: id_pot_rho_0 =-1 integer :: id_pot_rho_1 =-1 integer :: id_pot_rho_2 =-1 integer :: id_pot_rho_3 =-1 integer :: id_pot_rho_4 =-1 integer :: id_pot_rho_et =-1 integer :: id_pot_rho_nt =-1 integer :: id_pot_rho_wt =-1 integer :: id_rho_average =-1 integer :: id_mass_level =-1 integer :: id_volume_level =-1 integer :: id_drhodx_zt =-1 integer :: id_drhody_zt =-1 integer :: id_drhodz_zt =-1 integer :: id_drhodz_wt =-1 integer :: id_drhodz_diag =-1 integer :: id_buoyfreq2_zt =-1 integer :: id_buoyfreq2_wt =-1 integer :: id_cabbeling =-1 integer :: id_thermobaricity =-1 integer :: id_int_rhodz =-1 integer :: id_rhoave =-1 integer :: id_pbot_adjust =-1 integer :: id_eta_adjust =-1 integer :: id_eta_adjust_approx =-1 integer :: id_eta_adjust_cstvolume=-1 integer :: id_grad_nrho =-1 integer :: id_grad_lrpotrho =-1 integer :: id_grad_nrho_lrpotrho =-1 integer :: id_watermass_factor =-1 integer :: id_stratification_factor =-1 integer :: id_stratification_axis =-1 integer :: id_smax_dianeutral =-1 !for restart integer :: id_restart_rho = 0 integer :: id_restart_rho_s = 0 type(restart_file_type), save :: Den_restart #include type(ocean_domain_type), pointer :: Dom =>NULL() type(ocean_grid_type), pointer :: Grd =>NULL() character(len=128) :: version=& '$Id: ocean_density.F90,v 20.0 2013/12/14 00:10:40 fms Exp $' character (len=128) :: tagname = & '$Name: tikal $' public ocean_density_init public ocean_density_end public ocean_density_diag public update_ocean_density public update_ocean_density_salinity public density public potential_density public neutral_density public density_sfc public density_derivs public density_delta_z public density_delta_sfc public ocean_density_restart public calc_cabbeling_thermobaricity public buoyfreq2 public compute_buoyfreq private compute_drhodxy private ocean_density_chksum private cabbeling_thermobaricity private compute_diagnostic_factors private compute_density_diagnostics private density_coeffs_init private density_diagnostics_init ! interfaces for density interface density module procedure density_field module procedure density_level module procedure density_line module procedure density_point end interface ! interfaces for neutral density interface neutral_density module procedure neutral_density_field module procedure neutral_density_point end interface ! interfaces for density derivatives interface density_derivs module procedure density_derivs_field module procedure density_derivs_level module procedure density_derivs_point end interface integer :: index_temp =-1 integer :: index_salt =-1 integer :: index_Fdelta =-1 integer :: num_prog_tracers integer :: salt_variable integer :: temp_variable ! indices for diagnostic tracer density elements integer :: ind_rho = -1 integer :: ind_neutralrho = -1 integer :: ind_potrho_0 = -1 integer :: ind_potrho_1 = -1 integer :: ind_potrho_2 = -1 integer :: ind_potrho_3 = -1 integer :: ind_potrho_4 = -1 ! nml settings ! for linear EOS logical :: eos_linear=.false. real :: alpha_linear_eos=0.255 real :: beta_linear_eos =0.0 ! for teos10 or preteos10 eos logical :: eos_preteos10 = .false. logical :: eos_teos10 = .false. ! some test settings for density real :: s_test=20.0 real :: t_test=20.0 real :: p_test=1000.0 real :: jmfwg_rho real :: jmfwg_alpha real :: jmfwg_beta real :: mbfj_rho real :: mbfj_alpha real :: mbfj_beta ! some test settings for neutral density real :: sn_test=35.0 real :: tn_test=20.0 real :: mj_neutralrho real :: mb_neutralrho ! for diagnostic partitioning of vertical according to ! potential density, neutral density, or theta classes integer :: layer_nk = 80 ! # of classes used to compute diagnostics associated with ! potential density, potential temperature, or conservative ! temperature classes. ! for diagnostic partitioning of vertical according to potential density classes real :: potrho_press = 2000.0 ! sea pressure (dbars) for potential density computed for diagnostics real :: potrho_min = 1028.0 ! (kg/m^3) real :: potrho_max = 1038.0 ! (kg/m^3) ! for diagnostic partitioning of vertical according to ! rational polynomial approximation to neutral density logical :: neutral_density_omega = .false. logical :: neutral_density_potrho = .true. logical :: neutral_density_theta = .false. real :: neutralrho_min = 1020.0 ! (kg/m^3) real :: neutralrho_max = 1030.0 ! (kg/m^3) real :: nrho_press = 2000.0 logical :: diff_nrho_press = .false. ! for diagnostic partitioning of vertical according ! to potential temperature or conservative temperature classes real :: theta_min = -2.0 ! (degrees C) real :: theta_max = 30.0 ! (degrees C) ! standard atmospheric pressure (dbar). ! Should be set to 0.0 if assume zero pressure for the ! overlying atmosphere. But if have a realistic atmospheric ! pressure loading, then set press_standard=10.1325. real :: press_standard = 0.0 ! for smoothing the vertical density gradient ! when diagnosing the buoyancy frequency. logical :: buoyfreq_smooth_vert=.true. integer :: num_121_passes =1 real :: epsln_drhodz =1.e-10 ! for enforcing negative drhodz_diag logical :: drhodz_diag_stable = .true. real :: epsln_drhodz_diag = 1.e-10 ! for setting grad_nrho_lrpotrho=1.0 logical :: grad_nrho_lrpotrho_compute=.false. real :: grad_nrho_lrpotrho_max=10.0 real :: grad_nrho_lrpotrho_min=1.0 ! for stratification factor used for diagnostics logical :: smooth_stratification_factor=.false. logical :: smax_min_in_column=.false. real :: smax_diag=-1.0 real, dimension(:,:,:), allocatable :: smax real, dimension(:,:), allocatable :: dhorz_r ! to compute diagnostic factors for ! water mass transformation diagnostics logical :: update_diagnostic_factors=.false. ! for case when use domain masking, logical :: mask_domain_restart = .false. logical :: debug_this_module = .false. logical :: rho0_density = .false. logical :: density_equal_potrho = .false. logical :: module_is_initialized = .false. logical :: write_a_restart = .true. logical :: do_bitwise_exact_sum = .false. namelist /ocean_density_nml/ s_test, t_test, p_test, press_standard, & sn_test, tn_test, & eos_linear, alpha_linear_eos, beta_linear_eos, & eos_preteos10, eos_teos10, & potrho_press, potrho_min, potrho_max, & neutralrho_min, neutralrho_max, & layer_nk, theta_min, theta_max, & debug_this_module, write_a_restart, & rho0_density, density_equal_potrho, & buoyfreq_smooth_vert, num_121_passes, epsln_drhodz, & mask_domain_restart, do_bitwise_exact_sum, & drhodz_diag_stable, epsln_drhodz_diag, & grad_nrho_lrpotrho_compute, grad_nrho_lrpotrho_max, & grad_nrho_lrpotrho_min, & neutral_density_omega, neutral_density_potrho, & neutral_density_theta, & smooth_stratification_factor, update_diagnostic_factors, & smax_diag, smax_min_in_column contains !####################################################################### ! ! ! ! Initialize the density module ! ! subroutine ocean_density_init (Grid, Domain, Time, Time_steps, Thickness, T_prog, T_diag, & Ocean_options, Dens, ver_coordinate, use_blobs, debug) type(ocean_grid_type), target, intent(in) :: Grid type(ocean_domain_type), target, intent(in) :: Domain type(ocean_time_type), intent(in) :: Time type(ocean_time_steps_type), intent(in) :: Time_steps type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), target, intent(in) :: T_prog(:) type(ocean_diag_tracer_type), intent(inout) :: T_diag(:) type(ocean_options_type), intent(inout) :: Ocean_options type(ocean_density_type), target, intent(inout) :: Dens integer, intent(in) :: ver_coordinate logical, intent(in) :: use_blobs logical, intent(in), optional :: debug integer :: ioun, io_status, ierr integer :: i, j, k, n integer :: tau, taup1, taum1 integer :: num_eos real :: rho0_nk(nk) real, allocatable, dimension(:) :: rho_average real, allocatable, dimension(:) :: mass_level real, allocatable, dimension(:) :: volume_level integer :: id_grid_xt, id_grid_yt integer :: id_grid_xu, id_grid_yu integer :: id_potrho_bounds, id_potrho_axis integer :: id_neutralrho_bounds, id_neutralrho_axis integer :: id_theta_bounds, id_theta_axis integer :: id_restart(5) integer :: neutral_rho_method=0 real, allocatable, dimension(:) :: potrho_bounds real, allocatable, dimension(:) :: neutralrho_bounds real, allocatable, dimension(:) :: theta_bounds ! salinity and temperature pointers to aid readability real, dimension(:,:,:), pointer :: salinity => NULL() real, dimension(:,:,:), pointer :: temperature => NULL() real :: potrho_interval real :: neutralrho_interval real :: theta_interval character*128 filename integer :: stdoutunit,stdlogunit stdoutunit=stdout();stdlogunit=stdlog() if ( module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (ocean_density_init): module already initialized.') endif module_is_initialized = .TRUE. tau = Time%tau taup1 = Time%taup1 taum1 = Time%taum1 vert_coordinate = ver_coordinate num_prog_tracers = size(T_prog(:)) do n=1,num_prog_tracers if (T_prog(n)%name == 'temp') index_temp = n if (T_prog(n)%name == 'salt') index_salt = n if (T_prog(n)%name == 'fdelta') index_Fdelta = n if (T_prog(n)%longname == 'Conservative temperature') temp_variable = CONSERVATIVE_TEMP if (T_prog(n)%longname == 'Potential temperature') temp_variable = POTENTIAL_TEMP if (T_prog(n)%longname == 'Practical Salinity') salt_variable = PRACTICAL_SALT if (T_prog(n)%longname == 'Preformed Salinity') salt_variable = PREFORMED_SALT enddo call write_version_number(version, tagname) ! provide for namelist override of defaults #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=ocean_density_nml, iostat=io_status) ierr = check_nml_error(io_status,'ocean_density_nml') #else ioun = open_namelist_file() read (ioun,ocean_density_nml,IOSTAT=io_status) ierr = check_nml_error(io_status, 'ocean_density_nml') call close_file(ioun) #endif write (stdoutunit,'(/)') write (stdoutunit,ocean_density_nml) write (stdlogunit,ocean_density_nml) num_eos = 0 if(eos_linear) then num_eos = num_eos+1 elseif(eos_preteos10) then num_eos = num_eos+1 elseif(eos_teos10) then num_eos = num_eos+1 endif if(num_eos == 0) then call mpp_error(FATAL, & '==>Error in ocean_density_mod: No EOS chosen. Set one to true: eos_linear, eos_preteos10, eos_teso10.') endif if(num_eos > 1) then call mpp_error(FATAL, & '==>Error in ocean_density_mod: More than one EOS chosen. Set just one true: eos_linear, eos_preteos10, eos_teso10.') endif Dens%use_teos10 = eos_teos10 if (PRESENT(debug) .and. .not. debug_this_module) then debug_this_module = debug endif if(debug_this_module) then write(stdoutunit,'(a)') '==>Note: running ocean_density with debug_this_module=.true.' endif if(.not. write_a_restart) then write(stdoutunit,'(a)') '==>Note: running ocean_density with write_a_restart=.false.' write(stdoutunit,'(a)') ' Will NOT write restart file, and so cannot restart the run.' endif if(do_bitwise_exact_sum) then global_sum_flag = BITWISE_EXACT_SUM else global_sum_flag = NON_BITWISE_EXACT_SUM endif if(eos_linear) then write (stdoutunit,'(/1x,a)') & ' ==> Note: USING linear EOS designed for idealized simulations.' write (stdoutunit,'(7x,a)') & 'It is a linear function of potential temperature and salinity.' write (stdoutunit,'(7x,a)') & 'There is no pressure dependence.' Ocean_options%equation_of_state = 'Linear equation of state used for density.' elseif(eos_preteos10) then write (stdoutunit,'(/1x,a)') & ' ==> Note: USING preTEOS10 EOS, as relevant for realistic ocean climate simulations.' write (stdoutunit,'(1x,a,f12.6,a)') & ' Subtracting standard atmosphere of ',press_standard,' dbar for EOS calculation.' Ocean_options%equation_of_state = 'preTEOS10 equation of state used for density.' elseif(eos_teos10) then write (stdoutunit,'(/1x,a)') & ' ==> Note: USING TEOS10 EOS, as relevant for realistic ocean climate simulations.' write (stdoutunit,'(1x,a,f12.6,a)') & ' Subtracting standard atmosphere of ',press_standard,' dbar for EOS calculation.' Ocean_options%equation_of_state = 'TEOS10 equation of state used for density.' endif if(temp_variable==CONSERVATIVE_TEMP) then write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing EOS assuming prognostic temp = conservative temperature.' elseif(temp_variable==POTENTIAL_TEMP) then write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing EOS assuming prognostic temp = potential temperature.' else call mpp_error(FATAL, & '==>Error in ocean_density_mod: model temperature variable remains unspecified') endif if(salt_variable==PRACTICAL_SALT) then if(eos_teos10) then call mpp_error(FATAL, & '==>Error in ocean_density_mod: eos_teos10 requires preformed salinity as prognostic salinity.') endif write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing EOS assuming prognostic salinity = practical salinity.' elseif(salt_variable==PREFORMED_SALT) then if(eos_linear .or. eos_preteos10) then call mpp_error(FATAL, & '==>Error in ocean_density_mod: preformed salinity as prognostic salinity requires eos_teos10.') endif write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing TEOS10 EOS assuming prognostic salinity = preformed salinity.' else call mpp_error(FATAL, & '==>Error in ocean_density_mod: model salinity variable remains unspecified') endif write(stdoutunit,*) & '==>Note: The Boussinesq rho0 density has a value of (kg/m^3) ',rho0 if(rho0_density) then write(stdoutunit,*) & '==>Warning: rho0_density=.true, so rho=rho0 everywhere. Is this really what you wish to do?' endif if(neutral_density_omega) then write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing diagnostic neutral_rho according to Klocker/McDougall (not yet coded).' neutral_rho_method = neutral_rho_method + 1 call mpp_error(WARNING, & '==>Warning in ocean_density_mod: neutral_density_omega not coded. Will use potential density 2000dbar.') endif if(neutral_density_potrho) then write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing diagnostic neutral_rho as potential density referenced to pressure potrho_press.' neutral_rho_method = neutral_rho_method + 1 endif if(neutral_density_theta) then write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing diagnostic neutral_rho as conservative temperature.' neutral_rho_method = neutral_rho_method + 1 endif if(neutral_rho_method==0) then call mpp_error(FATAL, & '==>Error in ocean_density_mod: set one of the options "neutral_density_omega" or "neutral_density_potrho" or "neutral_density_theta" true.') endif if(neutral_rho_method > 1) then call mpp_error(FATAL, & '==>Error in ocean_density_mod: set only **one** of the options "neutral_density_omega" or "neutral_density_potrho" or "neutral_density_theta" true.') endif if(drhodz_diag_stable) then write (stdoutunit,'(/1x,a)') & ' ==> Note: Enforcing drhodz_diag < 0, so to use stable stratification for certain diagnostic purposes.' endif if(update_diagnostic_factors) then write (stdoutunit,'(/1x,a)') & ' ==> Note: Computing diagnostic factors in ocean_density.F90 for watermass diagnostics.' else write (stdoutunit,'(/1x,a)') & ' ==> Note: Diagnostic factors are NOT computed. So if enable watermass diagnostics, they will be corrupted.' endif #ifndef MOM_STATIC_ARRAYS call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) nk = Grid%nk allocate(Dens%rho(isd:ied,jsd:jed,nk,3)) allocate(Dens%rho_fresh(isd:ied,jsd:jed,nk)) allocate(Dens%rho_salinity(isd:ied,jsd:jed,nk,3)) allocate(Dens%rho_dztr_tau(isd:ied,jsd:jed,nk)) allocate(Dens%potrho(isd:ied,jsd:jed,nk)) allocate(Dens%neutralrho(isd:ied,jsd:jed,nk)) allocate(Dens%pressure_at_depth(isd:ied,jsd:jed,nk)) allocate(Dens%drhodT(isd:ied,jsd:jed,nk)) allocate(Dens%dpotrhodT(isd:ied,jsd:jed,nk)) allocate(Dens%dpotrhodS(isd:ied,jsd:jed,nk)) allocate(Dens%dpotrhodP(isd:ied,jsd:jed,nk)) allocate(Dens%drhodS(isd:ied,jsd:jed,nk)) allocate(Dens%drhodP(isd:ied,jsd:jed,nk)) allocate(Dens%drhodz_wt(isd:ied,jsd:jed,nk)) allocate(Dens%drhodz_zt(isd:ied,jsd:jed,nk)) allocate(Dens%drhodx_zt(isd:ied,jsd:jed,nk)) allocate(Dens%drhody_zt(isd:ied,jsd:jed,nk)) allocate(Dens%drhodz_diag(isd:ied,jsd:jed,nk)) allocate(Dens%watermass_factor(isd:ied,jsd:jed,nk)) allocate(Dens%stratification_factor(isd:ied,jsd:jed,nk)) allocate(Dens%dTdz_zt(isd:ied,jsd:jed,nk)) allocate(Dens%dSdz_zt(isd:ied,jsd:jed,nk)) allocate(Dens%mld_subduction(isd:ied,jsd:jed)) if (use_blobs) allocate(Dens%rhoT(isd:ied,jsd:jed,nk)) #endif do k=1,nk Dens%rho_dztr_tau(:,:,k) = rho0r*Grid%tmask(:,:,k) Dens%rho_fresh(:,:,k) = 1000.0*Grid%tmask(:,:,k) Dens%potrho(:,:,k) = rho0*Grid%tmask(:,:,k) Dens%neutralrho(:,:,k) = rho0*Grid%tmask(:,:,k) Dens%pressure_at_depth(:,:,k) = rho0*grav*Thickness%depth_zt(:,:,k)*c2dbars enddo Dens%rho(:,:,:,1) = Grid%tmask(:,:,:)*rho0 Dens%rho(:,:,:,2) = Grid%tmask(:,:,:)*rho0 Dens%rho(:,:,:,3) = Grid%tmask(:,:,:)*rho0 Dens%rho_salinity(:,:,:,1) = Grid%tmask(:,:,:)*35 Dens%rho_salinity(:,:,:,2) = Grid%tmask(:,:,:)*35 Dens%rho_salinity(:,:,:,3) = Grid%tmask(:,:,:)*35 Dens%drhodT(:,:,:) = 0.0 Dens%dpotrhodT(:,:,:) = 0.0 Dens%dpotrhodS(:,:,:) = 0.0 Dens%dpotrhodP(:,:,:) = 0.0 Dens%drhodS(:,:,:) = 0.0 Dens%drhodP(:,:,:) = 0.0 Dens%drhodz_wt(:,:,:) = Grid%tmask(:,:,:)*epsln_drhodz Dens%drhodz_zt(:,:,:) = Grid%tmask(:,:,:)*epsln_drhodz Dens%drhodx_zt(:,:,:) = 0.0 Dens%drhody_zt(:,:,:) = 0.0 Dens%drhodz_diag(:,:,:) = Grid%tmask(:,:,:)*epsln_drhodz_diag Dens%dTdz_zt(:,:,:) = 0.0 Dens%dSdz_zt(:,:,:) = 0.0 Dens%watermass_factor(:,:,:) = Grid%tmask(:,:,:) Dens%mld_subduction(:,:) = 0.0 do k=1,nk Dens%stratification_factor(:,:,k) = Grid%tmask(:,:,k)*Grid%dzt(k)*Grid%dat(:,:) enddo if (use_blobs) Dens%rhoT(:,:,:) = 0.0 Dom => Domain Grd => Grid dtts = Time_steps%dtts ! for diagnostics area_total_r = mpp_global_sum(Dom%domain2d,Grd%dat(:,:)*Grd%tmask(:,:,1), global_sum_flag) if (global_sum_flag .eq. NON_BITWISE_EXACT_SUM) area_total_r = real(area_total_r,kind=FLOAT_KIND) area_total_r = 1.0/(epsln+area_total_r) ! for diagnostic calculation of dianeutral transport allocate(smax(isd:ied,jsd:jed,nk)) allocate(dhorz_r(isd:ied,jsd:jed)) dhorz_r(:,:) = 0.0 smax(:,:,:) = 0.0 do j=jsd,jed do i=isd,ied dhorz_r(i,j) = Grd%tmask(i,j,1) & *( Grd%dxt(i,j) + Grd%dyt(i,j) ) / (epsln + 2.0*Grd%dxt(i,j)*Grd%dyt(i,j) ) enddo enddo grav_r = 1.0/grav allocate(rhodz_tau(isd:ied,jsd:jed)) allocate(rhodz_taup1(isd:ied,jsd:jed)) rhodz_tau = 0.0 rhodz_taup1 = 0.0 ! send as 1-dim field, due to problems with diag manager id_rho0 = register_static_field('ocean_model','rho0', Grid%tracer_axes(3:3), & 'reference density for Boussinesq approximation','kg/m^3', & missing_value=missing_value, range=(/-1e0,1e10/), & standard_name='reference_sea_water_density_for_boussinesq_approximation') rho0_nk(:) = rho0 if (id_rho0 > 0) then used = send_data(id_rho0, rho0_nk(:), Time%model_time) endif allocate(denominator_r(isd:ied,jsd:jed,nk)) denominator_r(:,:,:) = 0.0 ! specify coefficients for the polynomical equations of state ! and perform some pointwise tests. call density_coeffs_init() filename = 'ocean_density.res.nc' id_restart_rho = register_restart_field(Den_restart, filename, 'rho', Dens%rho(:,:,:,taup1), & domain=Dom%domain2d ) id_restart_rho_s = register_restart_field(Den_restart, filename, 'rho_salinity', Dens%rho_salinity(:,:,:,taup1), & domain=Dom%domain2d ) id_restart(1) = register_restart_field(Den_restart, filename, 'pressure_at_depth', Dens%pressure_at_depth(:,:,:), & domain=Dom%domain2d ) id_restart(2) = register_restart_field(Den_restart, filename, 'denominator_r', denominator_r(:,:,:), & domain=Dom%domain2d ) id_restart(3) = register_restart_field(Den_restart, filename, 'drhodT', Dens%drhodT(:,:,:), & domain=Dom%domain2d ) id_restart(4) = register_restart_field(Den_restart, filename, 'drhodS', Dens%drhodS(:,:,:), & domain=Dom%domain2d ) id_restart(5) = register_restart_field(Den_restart, filename, 'drhodz_zt', Dens%drhodz_zt(:,:,:), & domain=Dom%domain2d ) ! pointers to relevant salinity and temperature to enhance readability salinity => Dens%rho_salinity(:,:,:,tau) temperature => T_prog(index_temp)%field(:,:,:,tau) filename = 'INPUT/ocean_density.res.nc' if (.NOT.file_exist(trim(filename)) )then write (stdoutunit,'(/a)') ' Initialising salinity for use in density calculation' call update_ocean_density_salinity(T_prog,taum1,Dens) call update_ocean_density_salinity(T_prog,tau,Dens) call update_ocean_density_salinity(T_prog,taup1,Dens) Dens%rho(:,:,:,taum1) = density(salinity, temperature, Dens%pressure_at_depth(:,:,:)) Dens%rho(:,:,:,tau) = Dens%rho(:,:,:,taum1) Dens%rho(:,:,:,taup1) = Dens%rho(:,:,:,taum1) if (use_blobs) Dens%rhoT(:,:,:) = Dens%rho(:,:,:,taup1) else write (stdoutunit,'(/a)') ' Reading restart for density from INPUT/ocean_density.res.nc' call restore_state( Den_restart, id_restart_rho ) call mpp_update_domains(Dens%rho(:,:,:,taup1), Dom%domain2d) call restore_state( Den_restart, id_restart(1) ) call mpp_update_domains(Dens%pressure_at_depth(:,:,:),Dom%domain2d) call restore_state( Den_restart, id_restart(2) ) ! initialize to epsln to eliminate NaNs when mask ! out processors and then change layout upon a restart if(mask_domain_restart) then where (denominator_r==0.) denominator_r=rho0r endif call mpp_update_domains(denominator_r(:,:,:),Dom%domain2d) ! determine whether to read density derivs at the initial time step of integration. ! early versions of MOM4p1 did not contain drhodT and drhodS in the restart file, so ! to allow older restart files to be used with newer releases of MOM4p1, we ! provide for the possibility that the restarts do not contain drhodT and drhodS. if(field_exist(filename, 'drhodT')) then call restore_state( Den_restart, id_restart(3) ) call restore_state( Den_restart, id_restart(4) ) call mpp_update_domains(Dens%drhodT(:,:,:),Dom%domain2d) call mpp_update_domains(Dens%drhodS(:,:,:),Dom%domain2d) else write(stdoutunit,'(a)') '==>Note: ocean_density_mod: did not read density derivatives from restart.' endif ! determine whether to read rho_salinity at the initial time step of ! integration. Early versions of MOM4p1 did not contain the rho_salinity ! in the restart file. This allows older restart files to be used with ! newer releases of MOM. if(field_exist(filename, 'rho_salinity')) then write (stdoutunit,'(/a)') ' Initializing salinity for use in density calculation from restart' call restore_state( Den_restart, id_restart_rho_s ) Dens%rho_salinity(:,:,:,taum1) = Dens%rho_salinity(:,:,:,taup1) Dens%rho_salinity(:,:,:,tau) = Dens%rho_salinity(:,:,:,taup1) call mpp_update_domains(Dens%rho_salinity(:,:,:,taum1),Dom%domain2d) call mpp_update_domains(Dens%rho_salinity(:,:,:,tau), Dom%domain2d) call mpp_update_domains(Dens%rho_salinity(:,:,:,taup1),Dom%domain2d) else write (stdoutunit,'(/a)') ' Initialising salinity for use in density calculation' call update_ocean_density_salinity(T_prog,taum1,Dens) call update_ocean_density_salinity(T_prog,tau,Dens) call update_ocean_density_salinity(T_prog,taup1,Dens) endif endif ! compute buoyancy frequency for diagnostic purposes ! note: Dens%rhoT and the buoyancy frequency is initialised in ocean_blob_init ! when use_blobs=.true. if (.not. use_blobs) then call compute_buoyfreq(Time, Thickness, salinity, temperature, Dens, use_blobs) endif ! over-write drhodz_zt from compute_buoyfreq with drhodz_zt from restart. ! Note that we determine whether to read drhodz_zt at the initial time ! step of integration. early versions of MOM4p1 did not contain drhodz_zt ! in the restart file, so to allow older restart files to be used with ! newer releases of MOM4p1, we provide for the possibility that the ! restarts do not contain drhodz_zt. if (file_exist(trim(filename)) )then if(field_exist(filename, 'drhodz_zt')) then call restore_state( Den_restart, id_restart(5) ) call mpp_update_domains(Dens%drhodz_zt(:,:,:),Dom%domain2d) else write(stdoutunit,'(a)') '==>Note: ocean_density_mod: did not read drhodz_zt from restart.' endif endif ! compute neutral density for diagnostic purposes Dens%neutralrho(:,:,:) = Grd%tmask(:,:,:)*neutral_density(salinity, temperature) ind_neutralrho = fm_get_index('/ocean_mod/diag_tracers/rho_neutral') if (ind_neutralrho .gt. 0) then T_diag(ind_neutralrho)%field(:,:,:) = Dens%neutralrho(:,:,:) endif ! compute some density related diagnostic factors ! initialize neutralrho_interval_r first neutralrho_interval = (neutralrho_max-neutralrho_min)/(epsln+layer_nk) neutralrho_interval_r = 1.0/neutralrho_interval call compute_diagnostic_factors(Time, Thickness, Dens, salinity, temperature) ! compute potential density for diagnostic purposes Dens%potrho(:,:,:) = Grd%tmask(:,:,:)*potential_density(salinity, temperature, potrho_press) ! compute various potential densities for diagnostic purposes ind_potrho_0 = fm_get_index('/ocean_mod/diag_tracers/rho_pot_0') if (ind_potrho_0 .gt. 0) then T_diag(ind_potrho_0)%field(:,:,:) = Grd%tmask(:,:,:)*potential_density(salinity, temperature, 0.0) endif ind_potrho_1 = fm_get_index('/ocean_mod/diag_tracers/rho_pot_1') if (ind_potrho_1 .gt. 0) then T_diag(ind_potrho_1)%field(:,:,:) = Grd%tmask(:,:,:)*potential_density(salinity, temperature, 1000.0) endif ind_potrho_2 = fm_get_index('/ocean_mod/diag_tracers/rho_pot_2') if (ind_potrho_2 .gt. 0) then T_diag(ind_potrho_2)%field(:,:,:) = Grd%tmask(:,:,:)*potential_density(salinity, temperature, 2000.0) endif ind_potrho_3 = fm_get_index('/ocean_mod/diag_tracers/rho_pot_3') if (ind_potrho_3 .gt. 0) then T_diag(ind_potrho_3)%field(:,:,:) = Grd%tmask(:,:,:)*potential_density(salinity, temperature, 3000.0) endif ind_potrho_4 = fm_get_index('/ocean_mod/diag_tracers/rho_pot_4') if (ind_potrho_4 .gt. 0) then T_diag(ind_potrho_4)%field(:,:,:) = Grd%tmask(:,:,:)*potential_density(salinity, temperature, 4000.0) endif if(rho0_density) then Dens%rho(:,:,:,taum1) = Grd%tmask(:,:,:)*rho0 Dens%rho(:,:,:,tau) = Grd%tmask(:,:,:)*rho0 Dens%rho(:,:,:,taup1) = Grd%tmask(:,:,:)*rho0 endif ind_rho = fm_get_index('/ocean_mod/diag_tracers/rho_in_situ') if (ind_rho .gt. 0) then T_diag(ind_rho)%field(:,:,:) = Dens%rho(:,:,:,tau) endif ! compute initial level mass, volume, and density allocate(rho_average(0:nk)) allocate(mass_level(0:nk)) allocate(volume_level(0:nk)) rho_average(:) = 0.0 mass_level(:) = 0.0 volume_level(:) = 0.0 wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 wrk3_2d(:,:) = 0.0 do k=1,nk wrk1_2d(:,:) = Grd%dat(:,:)*Grd%tmask(:,:,k) wrk2_2d(:,:) = wrk1_2d(:,:)*Thickness%dzt(:,:,k) wrk3_2d(:,:) = wrk1_2d(:,:)*Thickness%dzt(:,:,k)*Dens%rho(:,:,k,tau) volume_level(k) = mpp_global_sum(Dom%domain2d,wrk2_2d(:,:),global_sum_flag) mass_level(k) = mpp_global_sum(Dom%domain2d,wrk3_2d(:,:),global_sum_flag) if (global_sum_flag == NON_BITWISE_EXACT_SUM) then volume_level(k) = real(volume_level(k),kind=FLOAT_KIND) mass_level(k) = real(mass_level(k),kind=FLOAT_KIND) endif if(volume_level(k) > 0.0) rho_average(k) = mass_level(k)/volume_level(k) enddo do k=1,nk volume_level(0) = volume_level(0) + volume_level(k) mass_level(0) = mass_level(0) + mass_level(k) enddo if(volume_level(0) > 0.0) rho_average(0) = mass_level(0)/volume_level(0) write(stdoutunit,'(/a,e24.12)') & ' ==>Note: From ocean_density_mod: Boussinesq reference density rho0(kg/m3) = ',rho0 write(stdoutunit,'(a,e24.12)') & ' Initial rho_average(kg/m3) = ',rho_average(0) if(rho_average(0) /= rho0 .and. Time%init) then write(stdoutunit,'(a)') & ' Since rho0 .ne. rho_average, consider changing rho0 in' write(stdoutunit,'(a/)') & ' ocean_parameters.F90 to be equal to rho_average for better accuracy.' endif id_rho_average = register_static_field('ocean_model','rho_average', Grid%tracer_axes(3:3), & 'level average density at initial step of segment','kg/m^3', & missing_value=missing_value, range=(/-1e0,1e10/)) id_mass_level = register_static_field('ocean_model','mass_level', Grid%tracer_axes(3:3), & 'mass of levels at initial step of segment','kg', & missing_value=missing_value, range=(/-1e0,1e10/)) id_volume_level= register_static_field('ocean_model','volume_level', Grid%tracer_axes(3:3), & 'volume of levels at initial step of segment','m^3', & missing_value=missing_value, range=(/-1e0,1e10/)) if (id_rho_average > 0) then used = send_data(id_rho_average, rho_average(1:nk), Time%model_time) endif if (id_mass_level > 0) then used = send_data(id_mass_level, mass_level(1:nk), Time%model_time) endif if (id_volume_level > 0) then used = send_data(id_volume_level, volume_level(1:nk), Time%model_time) endif ! define vertical axes according to potential density classes allocate ( Dens%potrho_ref(layer_nk)) allocate ( Dens%potrho_bounds(layer_nk+1)) allocate ( potrho_bounds(layer_nk+1)) potrho_bounds(1) = potrho_min potrho_interval = (potrho_max-potrho_min)/(epsln+layer_nk) do k=2,layer_nk+1 potrho_bounds(k)=potrho_bounds(k-1)+potrho_interval enddo do k=1,layer_nk Dens%potrho_ref(k)=potrho_bounds(k)+0.5*potrho_interval enddo do k=1,layer_nk+1 Dens%potrho_bounds(k) = potrho_bounds(k) enddo ! define vertical axes according to neutral density classes allocate ( Dens%neutralrho_ref(layer_nk)) allocate ( Dens%neutralrho_bounds(layer_nk+1)) allocate ( neutralrho_bounds(layer_nk+1)) neutralrho_bounds(1) = neutralrho_min do k=2,layer_nk+1 neutralrho_bounds(k)=neutralrho_bounds(k-1)+neutralrho_interval enddo do k=1,layer_nk Dens%neutralrho_ref(k)=neutralrho_bounds(k)+0.5*neutralrho_interval enddo do k=1,layer_nk+1 Dens%neutralrho_bounds(k) = neutralrho_bounds(k) enddo ! define vertical axes according to ! potential temperature or conservative temperature classes allocate ( Dens%theta_ref(layer_nk)) allocate ( Dens%theta_bounds(layer_nk+1)) allocate ( theta_bounds(layer_nk+1)) theta_bounds(1) = theta_min theta_interval = (theta_max-theta_min)/(epsln+layer_nk) do k=2,layer_nk+1 theta_bounds(k) =theta_bounds(k-1) +theta_interval enddo do k=1,layer_nk Dens%theta_ref(k) =theta_bounds(k) +0.5*theta_interval enddo do k=1,layer_nk+1 Dens%theta_bounds(k) = theta_bounds(k) enddo ! assuming that if the Grd%name = 'ocean' then we will ! be using this diagnostic axis. Otherwise it is not used. ! This is done to prevent duplicate axis names for multiple ! grid objects. if (trim(Grd%name) == 'ocean') then id_grid_xt = diag_axis_init ('grid_xt_'//trim(Grd%name),Grd%grid_x_t,'degrees_E','x', & 'tcell longitude',set_name='ocean', Domain2=Dom%domain2d) id_grid_yt = diag_axis_init ('grid_yt_'//trim(Grd%name),Grd%grid_y_t,'degrees_N','y', & 'tcell latitude',set_name='ocean', Domain2=Dom%domain2d) id_grid_xu = diag_axis_init ('grid_xu_'//trim(Grd%name),Grd%grid_x_u,'degrees_E','x', & 'ucell longitude',set_name='ocean', Domain2=Dom%domain2d) id_grid_yu = diag_axis_init ('grid_yu_'//trim(Grd%name),Grd%grid_y_u,'degrees_N','y', & 'ucell latitude',set_name='ocean', Domain2=Dom%domain2d) id_potrho_bounds = diag_axis_init & ( 'potrho_edges',potrho_bounds, 'kg/m^3', 'z','potential density edges',& direction=-1, set_name='ocean') id_potrho_axis = diag_axis_init & ( 'potrho',Dens%potrho_ref, & 'kg/m^3', 'z','potential density',edges=id_potrho_bounds,& direction=-1,set_name='ocean') Dens%potrho_axes = (/ id_grid_xt, id_grid_yt, id_potrho_axis /) Dens%potrho_axes_flux_x = (/ id_grid_xu, id_grid_yt, id_potrho_axis /) Dens%potrho_axes_flux_y = (/ id_grid_xt, id_grid_yu, id_potrho_axis /) id_neutralrho_bounds = diag_axis_init & ( 'neutralrho_edges',neutralrho_bounds, 'kg/m^3', 'z','neutral density edges',& direction=-1, set_name='ocean') id_neutralrho_axis = diag_axis_init & ( 'neutral',Dens%neutralrho_ref, & 'kg/m^3', 'z','neutral density',edges=id_neutralrho_bounds,& direction=-1,set_name='ocean') Dens%neutralrho_axes = (/ id_grid_xt, id_grid_yt, id_neutralrho_axis /) Dens%neutralrho_axes_flux_x = (/ id_grid_xu, id_grid_yt, id_neutralrho_axis /) Dens%neutralrho_axes_flux_y = (/ id_grid_xt, id_grid_yu, id_neutralrho_axis /) id_theta_bounds = diag_axis_init & ( 'theta_edges',potrho_bounds, 'C', 'z','potential or conservative temperature edges',& direction=1, set_name='ocean') id_theta_axis = diag_axis_init & ( 'theta',Dens%theta_ref, 'C', 'z','potential or conservative temperature',& edges=id_theta_bounds, direction=1,set_name='ocean') Dens%theta_axes = (/ id_grid_xt, id_grid_yt, id_theta_axis /) Dens%theta_axes_flux_x = (/ id_grid_xu, id_grid_yt, id_theta_axis /) Dens%theta_axes_flux_y = (/ id_grid_xt, id_grid_yu, id_theta_axis /) endif ! register the diagnostic fields call density_diagnostics_init(Time) write(stdoutunit,'(/a)') ' From ocean_density_mod: density chksums from ocean_density_init' call write_timestamp(Time%model_time) call ocean_density_chksum(Time, Dens, use_blobs) ! clean up pointers used in this subroutine nullify(salinity) nullify(temperature) end subroutine ocean_density_init ! NAME="ocean_density_init" !####################################################################### ! ! ! ! Register the diagnostic fields. ! ! subroutine density_diagnostics_init(Time) type(ocean_time_type), intent(in) :: Time id_press = register_diag_field ('ocean_model', 'press', Grd%tracer_axes(1:3),& Time%model_time, 'absolute pressure', 'dbar', & missing_value=missing_value, range=(/-10.,1e6/)) id_rho = register_diag_field ('ocean_model', 'rho', Grd%tracer_axes(1:3), & Time%model_time, 'in situ density', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/)) id_rhoE = register_diag_field ('ocean_model', 'rhoE', Grd%tracer_axes(1:3), & Time%model_time, 'E system in situ density', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/)) id_pot_rho = register_diag_field ('ocean_model', 'pot_rho', Grd%tracer_axes(1:3), & Time%model_time, 'potential density', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/)) id_int_rhodz = register_diag_field ('ocean_model', 'int_rhodz', Grd%tracer_axes(1:2), & Time%model_time, 'vertical integral of density', 'm*(kg/m^3)', & missing_value=missing_value, range=(/-10.0,1e5/)) if(neutral_density_omega) then id_neutral_rho = register_diag_field ('ocean_model', 'neutral_rho', Grd%tracer_axes(1:3), & Time%model_time, 'polynomial estimate of neutral density', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/)) elseif(neutral_density_potrho) then id_neutral_rho = register_diag_field ('ocean_model', 'neutral_rho', Grd%tracer_axes(1:3), & Time%model_time, 'potential density estimate of neutral density', 'kg/m^3',& missing_value=missing_value, range=(/-10.0,1e5/)) elseif(neutral_density_theta) then id_neutral_rho = register_diag_field ('ocean_model', 'neutral_rho', Grd%tracer_axes(1:3), & Time%model_time, 'conservative temperature is neutral density', 'degC',& missing_value=missing_value, range=(/-10.0,1e5/)) endif id_eos_salinity = register_diag_field ('ocean_model', 'eos_salinity', Grd%tracer_axes(1:3),& Time%model_time, 'absolute salinity as used for TEOS10 equation of state', 'g/kg',& missing_value=missing_value, range=(/0.0,100.0/)) id_pot_rho_0 = register_diag_field ('ocean_model', 'pot_rho_0', Grd%tracer_axes(1:3), & Time%model_time, 'potential density referenced to 0 dbar', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/), & standard_name='sea_water_potential_density') id_pot_rho_1 = register_diag_field ('ocean_model', 'pot_rho_1', Grd%tracer_axes(1:3), & Time%model_time, 'potential density referenced to 1000 dbar', 'kg/m^3',& missing_value=missing_value, range=(/-10.0,1e5/)) id_pot_rho_2 = register_diag_field ('ocean_model', 'pot_rho_2', Grd%tracer_axes(1:3), & Time%model_time, 'potential density referenced to 2000 dbar', 'kg/m^3',& missing_value=missing_value, range=(/-10.0,1e5/)) id_pot_rho_3 = register_diag_field ('ocean_model', 'pot_rho_3', Grd%tracer_axes(1:3), & Time%model_time, 'potential density referenced to 3000 dbar', 'kg/m^3',& missing_value=missing_value, range=(/-10.0,1e5/)) id_pot_rho_4 = register_diag_field ('ocean_model', 'pot_rho_4', Grd%tracer_axes(1:3), & Time%model_time, 'potential density referenced to 4000 dbar', 'kg/m^3',& missing_value=missing_value, range=(/-10.0,1e5/)) id_pot_rho_et = register_diag_field ('ocean_model', 'pot_rho_et', Grd%tracer_axes_flux_x(1:3), & Time%model_time, 'potential density at east face of tracer cell', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/)) id_pot_rho_nt = register_diag_field ('ocean_model', 'pot_rho_nt', Grd%tracer_axes_flux_y(1:3), & Time%model_time, 'potential density at north face of tracer cell', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/)) id_pot_rho_wt = register_diag_field ('ocean_model', 'pot_rho_wt', Grd%tracer_axes_wt(1:3), & Time%model_time, 'potential density at wt points', 'kg/m^3', & missing_value=missing_value, range=(/-10.0,1e5/)) id_thermal_expand = register_diag_field ('ocean_model', 'thermal_expand', Grd%tracer_axes(1:3), & Time%model_time, 'thermal expansion coefficient -(1/rho)d(rho)/d(theta)', '1/C',& missing_value=-1e10, range=(/-1e10,1e10/)) id_haline_contract= register_diag_field ('ocean_model', 'haline_contract', Grd%tracer_axes(1:3), & Time%model_time, 'haline contraction coefficient (1/rho)d(rho)/d(salinity)', '1/(g/kg)',& missing_value=-1e10, range=(/-1e10,1e10/)) id_drhodtheta = register_diag_field ('ocean_model', 'drhodtheta', Grd%tracer_axes(1:3), & Time%model_time, 'd(rho)/d(theta)', 'kg/m^3/C', & missing_value=-1e10, range=(/-1e10,1e10/)) id_drhodsalt = register_diag_field ('ocean_model', 'drhodsalinity', Grd%tracer_axes(1:3), & Time%model_time, 'd(rho)/d(salinity)', 'kg/m^3/psu', & missing_value=missing_value, range=(/-1e9,1e9/)) id_dpotrhodtheta = register_diag_field ('ocean_model', 'dpotrhodtheta', Grd%tracer_axes(1:3), & Time%model_time, 'd(potrho)/d(theta)', 'kg/m^3/C', & missing_value=-1e10, range=(/-1e10,1e10/)) id_dpotrhodsalt = register_diag_field ('ocean_model', 'dpotrhodsalt', Grd%tracer_axes(1:3), & Time%model_time, 'd(potrho)/d(salinity)', 'kg/m^3/psu', & missing_value=missing_value, range=(/-1e9,1e9/)) id_dpotrhodpress = register_diag_field ('ocean_model', 'dpotrhodpress', Grd%tracer_axes(1:3), & Time%model_time, 'd(potrho)/d(press)', 'kg/m^3/psu', & missing_value=missing_value, range=(/-1e9,1e9/)) id_drhodpress = register_diag_field ('ocean_model', 'drhodpress', Grd%tracer_axes(1:3), & Time%model_time, 'd(rho)/d(press)', '(kg/m^3)/dbar', & missing_value=missing_value, range=(/-1e9,1e9/)) id_sound_speed2 = register_diag_field ('ocean_model', 'sound_speed2', Grd%tracer_axes(1:3),& Time%model_time, 'squared sound speed = 1/[d(rho)/d(press)]', '(m/s)^2', & missing_value=missing_value, range=(/-1e9,1e9/)) id_cabbeling = register_diag_field ('ocean_model', 'cabbeling', Grd%tracer_axes(1:3), & Time%model_time, 'cabbeling parameter', '(1/degC)^2', & missing_value=-1e10, range=(/-1e10,1e10/)) id_thermobaricity = register_diag_field ('ocean_model', 'thermobaricity', Grd%tracer_axes(1:3), & Time%model_time, 'thermobaricity parameter', '1/(dbar*degC)', & missing_value=-1e10, range=(/-1e10,1e10/)) id_buoyfreq2_zt = register_diag_field ('ocean_model','buoyfreq2_zt', & Grd%tracer_axes(1:3), Time%model_time, & 'Squared buoyancy frequency at T-point', '1/s^2', & missing_value=missing_value, range=(/-1e6,1e6/)) id_buoyfreq2_wt = register_diag_field ('ocean_model','buoyfreq2_wt', & Grd%tracer_axes_wt(1:3), Time%model_time, & 'Squared buoyancy frequency at T-cell bottom', '1/s^2', & missing_value=missing_value, range=(/-1e6,1e6/)) id_drhodz_zt = register_diag_field ('ocean_model','drhodz_zt',& Grd%tracer_axes(1:3), Time%model_time, & 'd(neutral rho)/dz at T-point', 'kg/m^4', & missing_value=missing_value, range=(/-1e6,1e6/)) id_drhodx_zt = register_diag_field ('ocean_model','drhodx_zt', & Grd%tracer_axes(1:3), Time%model_time, & 'd(neutral rho)/dx at nominal T-point', 'kg/m^4',& missing_value=missing_value, range=(/-1e6,1e6/)) id_drhody_zt = register_diag_field ('ocean_model','drhody_zt', & Grd%tracer_axes(1:3), Time%model_time, & 'd(neutral rho)/dy at nominal T-point', 'kg/m^4',& missing_value=missing_value, range=(/-1e6,1e6/)) id_drhodz_diag = register_diag_field ('ocean_model','drhodz_diag', & Grd%tracer_axes(1:3), Time%model_time, & 'regularized d(neutral rho)/dz at T-point', 'kg/m^4',& missing_value=missing_value, range=(/-1e6,1e6/)) id_drhodz_wt = register_diag_field ('ocean_model','drhodz_wt', & Grd%tracer_axes_wt(1:3), Time%model_time, & 'd(neutral rho)/dz at W-point', 'kg/m^4', & missing_value=missing_value, range=(/-1e6,1e6/)) id_rhoave = register_diag_field('ocean_model','rhoave', Time%model_time, & 'global mean ocean in-situ density from ocean_density_mod',& 'kg/m^3', missing_value=missing_value, range=(/-1e10,1e10/)) id_pbot_adjust = register_diag_field('ocean_model','pbot_adjust',Time%model_time, & 'pbot adjustment to counteract spurious mass source in Boussinesq fluid',& 'dbar', missing_value=missing_value, range=(/-1e10,1e10/)) id_eta_adjust = register_diag_field('ocean_model','eta_adjust',Time%model_time, & 'global eta adjustment to include steric effect in Boussinesq fluid',& 'm', missing_value=missing_value, range=(/-100.0,100.0/)) id_eta_adjust_approx = register_diag_field('ocean_model','eta_adjust_approx',Time%model_time, & 'approximate global eta adjustment to include steric effect in Boussinesq fluid',& 'm', missing_value=missing_value, range=(/-100.0,100.0/)) id_eta_adjust_cstvolume = register_diag_field('ocean_model','eta_adjust_cstvolume',Time%model_time, & 'global eta adjustment (assuming constant volume) to include steric effect in Boussinesq fluid',& 'm', missing_value=missing_value, range=(/-100.0,100.0/)) id_grad_nrho = register_diag_field ('ocean_model','grad_nrho', & Grd%tracer_axes(1:3), Time%model_time, & '|grad neutral rho|', 'dimensionless', & missing_value=missing_value, range=(/-1e1,1e10/)) id_grad_lrpotrho = register_diag_field ('ocean_model','grad_lrpotrho', & Grd%tracer_axes(1:3), Time%model_time, & '|grad local ref potential rho|', 'dimensionless', & missing_value=missing_value, range=(/-1e1,1e10/)) id_grad_nrho_lrpotrho = register_diag_field ('ocean_model','grad_nrho_lrpotrho', & Grd%tracer_axes(1:3), Time%model_time, & '|grad neutral rho| / |grad local ref potential rho|', 'dimensionless',& missing_value=missing_value, range=(/-1e1,1e10/)) id_watermass_factor = register_diag_field ('ocean_model','watermass_factor', & Grd%tracer_axes(1:3), Time%model_time, & '(|grad neutral rho|/|grad local ref potential rho|)/(delta(gamma))', '(m^3/kg)',& missing_value=missing_value, range=(/-1e1,1e15/)) id_stratification_factor = register_diag_field ('ocean_model','stratification_factor', & Grd%tracer_axes(1:3), Time%model_time, 'abs ( rho*Area(h)/gamma_h )', 'm^3', & missing_value=missing_value, range=(/-1e1,1e22/)) id_stratification_axis = register_diag_field ('ocean_model','stratification_axis',Grd%tracer_axes(1:3), & Time%model_time, 'direction of strongest stratification', 'dimensionless', & missing_value=missing_value, range=(/-1e2,1e2/)) id_smax_dianeutral = register_diag_field ('ocean_model','smax_dianeutral',Grd%tracer_axes(1:3), & Time%model_time, 'slope factor used for computing dianeutral transport diagnostic',& 'dimensionless', missing_value=missing_value, range=(/-1e2,1e2/)) return end subroutine density_diagnostics_init ! NAME="density_diagnostics_init" !####################################################################### ! ! ! ! Initialize the EOS coefficients, and write some test values. ! ! subroutine density_coeffs_init() integer :: stdoutunit real :: rho_neutralrho real :: neutralrho_test, rho_test, alpha_test, beta_test, speed2_test real :: drho_dtheta_test, drho_dsal_test, drho_dpress_test real :: diff_test stdoutunit=stdout() ! for the TESO10 EOS mbfj_rho = 1.017775176234136d+3 mbfj_alpha = 2.435473441547041d-4 mbfj_beta = 7.284367916939847d-4 mb_neutralrho=1033.093610463980 v01 = 9.998420897506056d+2 v02 = 2.839940833161907 v03 = -3.147759265588511d-2 v04 = 1.181805545074306d-3 v05 = -6.698001071123802 v06 = -2.986498947203215d-2 v07 = 2.327859407479162d-4 v08 = -3.988822378968490d-2 v09 = 5.095422573880500d-4 v10 = -1.426984671633621d-5 v11 = 1.645039373682922d-7 v12 = -2.233269627352527d-2 v13 = -3.436090079851880d-4 v14 = 3.726050720345733d-6 v15 = -1.806789763745328d-4 v16 = 6.876837219536232d-7 v17 = -3.087032500374211d-7 v18 = -1.988366587925593d-8 v19 = -1.061519070296458d-11 v20 = 1.550932729220080d-10 v21 = 1.0 v22 = 2.775927747785646d-3 v23 = -2.349607444135925d-5 v24 = 1.119513357486743d-6 v25 = 6.743689325042773d-10 v26 = -7.521448093615448d-3 v27 = -2.764306979894411d-5 v28 = 1.262937315098546d-7 v29 = 9.527875081696435d-10 v30 = -1.811147201949891d-11 v31 = -3.303308871386421d-5 v32 = 3.801564588876298d-7 v33 = -7.672876869259043d-9 v34 = -4.634182341116144d-11 v35 = 2.681097235569143d-12 v36 = 5.419326551148740d-6 v37 = -2.742185394906099d-5 v38 = -3.212746477974189d-7 v39 = 3.191413910561627d-9 v40 = -1.931012931541776d-12 v41 = -1.105097577149576d-7 v42 = 6.211426728363857d-10 v43 = -1.119011592875110d-10 v44 = -1.941660213148725d-11 v45 = -1.864826425365600d-14 v46 = 1.119522344879478d-14 v47 = -1.200507748551599d-15 v48 = 6.057902487546866d-17 ! Save some multiples two_v03 = 2.0*v03 two_v07 = 2.0*v07 two_v10 = 2.0*v10 two_v14 = 2.0*v14 two_v17 = 2.0*v17 two_v18 = 2.0*v18 two_v19 = 2.0*v19 two_v20 = 2.0*v20 two_v23 = 2.0*v23 two_v28 = 2.0*v28 two_v33 = 2.0*v33 two_v36 = 2.0*v36 two_v39 = 2.0*v39 two_v43 = 2.0*v43 two_v44 = 2.0*v44 two_v45 = 2.0*v45 two_v46 = 2.0*v46 three_v04 = 3.0*v04 three_v11 = 3.0*v11 three_v24 = 3.0*v24 three_v29 = 3.0*v29 three_v34 = 3.0*v34 three_v40 = 3.0*v40 three_v47 = 3.0*v47 three_v48 = 3.0*v48 four_v25 = 4.0*v25 four_v30 = 4.0*v30 four_v35 = 4.0*v35 onep5_v08 = 1.5*v08 onep5_v09 = 1.5*v09 onep5_v10 = 1.5*v10 onep5_v11 = 1.5*v11 onep5_v10 = 1.5*v10 onep5_v31 = 1.5*v31 onep5_v32 = 1.5*v32 onep5_v33 = 1.5*v33 onep5_v34 = 1.5*v34 onep5_v35 = 1.5*v35 ! Double derivs six_v04 = 6.0*v04 six_v11 = 6.0*v11 six_v24 = 6.0*v24 six_v29 = 6.0*v29 six_v34 = 6.0*v34 six_v40 = 6.0*v40 six_v47 = 6.0*v47 six_v48 = 6.0*v48 twelve_v25 = 12.0*v25 twelve_v30 = 12.0*v30 twelve_v35 = 12.0*v35 p75_v08 = 0.75*v08 p75_v09 = 0.75*v09 p75_v10 = 0.75*v10 p75_v11 = 0.75*v11 p75_v10 = 0.75*v10 p75_v31 = 0.75*v31 p75_v32 = 0.75*v32 p75_v33 = 0.75*v33 p75_v34 = 0.75*v34 p75_v35 = 0.75*v35 ! Cross derivs three_v10 = 3.0*v10 three_v33 = 3.0*v33 fourp5_v11 = 4.5*v11 fourp5_v34 = 4.5*v34 four_v19 = 4.0*v19 four_v45 = 4.0*v45 six_v35 = 6.0*v35 ! 25 coefficients in the preTEOS10 equation of state if(temp_variable==CONSERVATIVE_TEMP) then jmfwg_rho = 1017.842890411975d0 jmfwg_alpha = 2.436057013634649d-4 jmfwg_beta = 7.314818108935248d-4 a0 = 9.9983912878771446d+02 a1 = 7.0687133522652896d+00 a2 = -2.2746841916232965d-02 a3 = 5.6569114861400121d-04 a4 = 2.3849975952593345d+00 a5 = 3.1761924314867009d-04 a6 = 1.7459053010547962d-03 a7 = 1.2192536310173776d-02 a8 = 2.4643435731663949d-07 a9 = 4.0525405332794888d-06 a10 = -2.3890831309113187d-08 a11 = -5.9016182471196891d-12 b0 = 1.0000000000000000d+00 b1 = 7.0051665739672298d-03 b2 = -1.5040804107377016d-05 b3 = 5.3943915288426715d-07 b4 = 3.3811600427083414d-10 b5 = 1.5599507046153769d-03 b6 = -1.8137352466500517d-06 b7 = -3.3580158763335367d-10 b8 = 5.7149997597561099d-06 b9 = 7.8025873978107375d-10 b10 = 7.1038052872522844d-06 b11 = -2.1692301739460094d-17 b12 = -8.2564080016458560d-18 ! Coefficients for neutral density based on McDougall/Jackett (2005). ! To be replaced by Klocker/McDougall approach in near future. rho_neutralrho=1024.43863927763d0 a0n = 1.0022048243661291d+003 a1n = 2.0634684367767725d-001 a2n = 8.0483030880783291d-002 a3n = -3.6670094757260206d-004 a4n = -1.4602011474139313d-003 a5n = -2.5860953752447594d-003 a6n = -3.0498135030851449d-007 b0n = 1.0000000000000000d+000 b1n = 4.4946117492521496d-005 b2n = 7.9275128750339643d-005 b3n = -1.2358702241599250d-007 b4n = -4.1775515358142458d-009 b5n = -4.3024523119324234d-004 b6n = 6.3377762448794933d-006 b7n = -7.2640466666916413d-010 b8n = -5.1075068249838284d-005 b9n = -5.8104725917890170d-009 elseif(temp_variable==POTENTIAL_TEMP) then jmfwg_rho = 1017.728868019642d0 jmfwg_alpha = 2.525481286927133d-4 jmfwg_beta = 7.379638527217575d-4 a0 = 9.9984085444849347d+02 a1 = 7.3471625860981584d+00 a2 = -5.3211231792841769d-02 a3 = 3.6492439109814549d-04 a4 = 2.5880571023991390d+00 a5 = -6.7168282786692355d-03 a6 = 1.9203202055760151d-03 a7 = 1.1798263740430364d-02 a8 = 9.8920219266399117d-08 a9 = 4.6996642771754730d-06 a10 = -2.5862187075154352d-08 a11 = -3.2921414007960662d-12 b0 = 1.0000000000000000d+00 b1 = 7.2815210113327091d-03 b2 = -4.4787265461983921d-05 b3 = 3.3851002965802430d-07 b4 = 1.3651202389758572d-10 b5 = 1.7632126669040377d-03 b6 = -8.8066583251206474d-06 b7 = -1.8832689434804897d-10 b8 = 5.7463776745432097d-06 b9 = 1.4716275472242334d-09 b10 = 6.7103246285651894d-06 b11 = -2.4461698007024582d-17 b12 = -9.1534417604289062d-18 ! Coefficients for neutral density based on McDougall/Jackett (2005). ! To be replaced by Klocker/McDougall approach in near future. rho_neutralrho=1024.59416751197d0 a0n = 1.0023063688892480d+003 a1n = 2.2280832068441331d-001 a2n = 8.1157118782170051d-002 a3n = -4.3159255086706703d-004 a4n = -1.0304537539692924d-004 a5n = -3.1710675488863952d-003 a6n = -1.7052298331414675d-007 b0n = 1.0000000000000000d+000 b1n = 4.3907692647825900d-005 b2n = 7.8717799560577725d-005 b3n = -1.6212552470310961d-007 b4n = -2.3850178558212048d-009 b5n = -5.1268124398160734d-004 b6n = 6.0399864718597388d-006 b7n = -2.2744455733317707d-009 b8n = -3.6138532339703262d-005 b9n = -1.3409379420216683d-009 endif ! save some multiples of the coefficients two_a2 = 2.0*a2 three_a3 = 3.0*a3 six_a3 = 6.0*a3 two_a6 = 2.0*a6 two_a8 = 2.0*a8 two_a10 = 2.0*a10 two_a11 = 2.0*a11 four_a11 = 4.0*a11 two_b2 = 2.0*b2 three_b3 = 3.0*b3 six_b3 = 6.0*b3 four_b4 = 4.0*b4 twelve_b4 = 12.0*b4 three_b7 = 3.0*b7 six_b7 = 6.0*b7 onep5_b8 = 1.5*b8 onep5_b9 = 1.5*b9 two_b9 = 2.0*b9 three_b9 = 3.0*b9 two_b11 = 2.0*b11 three_b11 = 3.0*b11 six_b11 = 6.0*b11 three_b12 = 3.0*b12 ! Test values for preTEOS10 EOS if(eos_preteos10) then write (stdoutunit,'(/,a)') 'preTEOS10 EQUATION OF STATE TEST VALUES' write (stdoutunit,'(a,f6.2,a,f6.2,a,f8.2)') & 's_test(psu) = ',s_test,', t_test(C) = ',t_test,', p_test(dbar) = ',p_test rho_test = density(s_test,t_test,p_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'rho (',s_test,',',t_test,',',p_test,') = ',rho_test,' kg/m^3' diff_test = rho_test-jmfwg_rho write (stdoutunit,' (a,e22.16,a)') 'diff from JMFWG = ',diff_test,' kg/m^3' call density_derivs(rho_test, s_test, t_test, p_test, & drho_dtheta_test, drho_dsal_test, drho_dpress_test) alpha_test = -drho_dtheta_test/(epsln+rho_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'alpha(',s_test,',',t_test,',',p_test,') = ',alpha_test,' 1/C' diff_test = alpha_test-jmfwg_alpha write (stdoutunit,' (a,e22.16,a)') 'diff from JMFWG = ',diff_test,' 1/C' beta_test = drho_dsal_test /(epsln+rho_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'beta (',s_test,',',t_test,',',p_test,') = ',beta_test,' 1/psu' speed2_test = 1.0/(epsln + drho_dpress_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'squared sound speed (',s_test,',',t_test,',',p_test,') = ',speed2_test,' (m/s)^2' diff_test = beta_test-jmfwg_beta write (stdoutunit,' (a,e22.16,a)') 'diff from JMFWG = ',diff_test,' 1/psu' write (stdoutunit,'(/,a)') 'NEUTRAL DENSITY EQUATION TEST VALUES' write (stdoutunit,'(a,f6.2,a,f6.2)') & 'sn_test(psu) = ',sn_test,', tn_test(C) = ',tn_test neutralrho_test = neutral_density(sn_test,tn_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,e22.16,a)') & 'rho (',sn_test,',',tn_test,') = ',neutralrho_test,' kg/m^3' diff_test = neutralrho_test-rho_neutralrho write (stdoutunit,' (a,e22.16,a)') 'diff from Klocker and McDougall test = ',diff_test,' kg/m^3' endif ! Test values for TEOS10 EOS if (eos_teos10) then write (stdoutunit,'(/,a)') 'TEOS10 EQUATION OF STATE TEST VALUES' write (stdoutunit,'(a,f6.2,a,f6.2,a,f8.2)') & 's_test(ppt) = ',s_test,', t_test(C) = ',t_test,', p_test(dbar) = ',p_test rho_test = density(s_test,t_test,p_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'rho (',s_test,',',t_test,',',p_test,') = ',rho_test,' kg/m^3' diff_test = rho_test-mbfj_rho write (stdoutunit,' (a,e22.16,a)') 'diff from MBFJ = ',diff_test,' kg/m^3' call density_derivs(rho_test, s_test, t_test, p_test, & drho_dtheta_test, drho_dsal_test, drho_dpress_test) alpha_test = -drho_dtheta_test/(epsln+rho_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'alpha(',s_test,',',t_test,',',p_test,') = ',alpha_test,' 1/C' diff_test = alpha_test-mbfj_alpha write (stdoutunit,' (a,e22.16,a)') 'diff from MBFJ = ',diff_test,' 1/C' beta_test = drho_dsal_test /(epsln+rho_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'beta (',s_test,',',t_test,',',p_test,') = ',beta_test,' 1/ppt' diff_test = beta_test-mbfj_beta write (stdoutunit,' (a,e22.16,a)') 'diff from MBFJ = ',diff_test,' 1/ppt' speed2_test = 1.0/(epsln + drho_dpress_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,f8.2,a,e22.16,a)') & 'squared sound speed (',s_test,',',t_test,',',p_test,') = ',speed2_test,' (m/s)^2' write (stdoutunit,'(/,a)') 'NEUTRAL DENSITY EQUATION TEST VALUES' write (stdoutunit,'(a,f6.2,a,f6.2)') & 'sn_test(ppt) = ',sn_test,', tn_test(C) = ',tn_test neutralrho_test = neutral_density(sn_test,tn_test) write (stdoutunit,' (a,f6.2,a,f6.2,a,e22.16,a)') & 'rho (',sn_test,',',tn_test,') = ',neutralrho_test,' kg/m^3' diff_test = neutralrho_test-mb_neutralrho write (stdoutunit,' (a,e22.16,a)') 'diff from McDougall and Barker = ',diff_test,' kg/m^3' endif return end subroutine density_coeffs_init ! NAME="density_coeffs_init" !####################################################################### ! ! ! ! Diagnostic ocean density fields: neutral density and potential density. ! Also send some diagnostics to diagnostic manager. ! ! subroutine ocean_density_diag(Time, Temp, Thickness, Dens, T_diag) type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), intent(in) :: Temp type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(inout) :: Dens type(ocean_diag_tracer_type), intent(inout) :: T_diag(:) integer :: tau tau = Time%tau if (ind_rho .gt. 0) then T_diag(ind_rho)%field(:,:,:) = Dens%rho(:,:,:,tau) endif ! compute neutral density for diagnostic purposes Dens%neutralrho(:,:,:) = Grd%tmask(:,:,:) & *neutral_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau)) if (ind_neutralrho .gt. 0) then T_diag(ind_neutralrho)%field(:,:,:) = Dens%neutralrho(:,:,:) endif ! compute potential density for diagnostic purposes Dens%potrho(:,:,:) = Grd%tmask(:,:,:) & *potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), potrho_press) ! compute some diagnostic factors call compute_diagnostic_factors(Time, Thickness, Dens, & Dens%rho_salinity(:,:,:,tau),& Temp%field(:,:,:,tau)) ! compute potential densities for diagnostic purposes if (ind_potrho_0 .gt. 0) then T_diag(ind_potrho_0)%field(:,:,:) = Grd%tmask(:,:,:) * & potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 0.0) endif if (ind_potrho_1 .gt. 0) then T_diag(ind_potrho_1)%field(:,:,:) = Grd%tmask(:,:,:) * & potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 1000.0) endif if (ind_potrho_2 .gt. 0) then T_diag(ind_potrho_2)%field(:,:,:) = Grd%tmask(:,:,:) * & potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 2000.0) endif if (ind_potrho_3 .gt. 0) then T_diag(ind_potrho_3)%field(:,:,:) = Grd%tmask(:,:,:) * & potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 3000.0) endif if (ind_potrho_4 .gt. 0) then T_diag(ind_potrho_4)%field(:,:,:) = Grd%tmask(:,:,:) * & potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 4000.0) endif call diagnose_3d(Time, Grd, id_neutral_rho, Dens%neutralrho(:,:,:)) call diagnose_3d(Time, Grd, id_pot_rho, Dens%potrho(:,:,:)) if (id_pot_rho_0 > 0) then if (ind_potrho_0 .gt. 0) then call diagnose_3d(Time, Grd, id_pot_rho_0, T_diag(ind_potrho_0)%field(:,:,:)) else wrk1(:,:,:) = potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 0.0) call diagnose_3d(Time, Grd, id_pot_rho_0, wrk1(:,:,:)) endif endif if (id_pot_rho_1 > 0) then if (ind_potrho_1 .gt. 0) then call diagnose_3d(Time, Grd, id_pot_rho_1, T_diag(ind_potrho_1)%field(:,:,:)) else wrk1(:,:,:) = potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 1000.0) call diagnose_3d(Time, Grd, id_pot_rho_1, wrk1(:,:,:)) endif endif if (id_pot_rho_2 > 0) then if (ind_potrho_2 .gt. 0) then call diagnose_3d(Time, Grd, id_pot_rho_2, T_diag(ind_potrho_2)%field(:,:,:)) else wrk1(:,:,:) = potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 2000.0) call diagnose_3d(Time, Grd, id_pot_rho_2, wrk1(:,:,:)) endif endif if (id_pot_rho_3 > 0) then if (ind_potrho_3 .gt. 0) then call diagnose_3d(Time, Grd, id_pot_rho_3, T_diag(ind_potrho_3)%field(:,:,:)) else wrk1(:,:,:) = potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 3000.0) call diagnose_3d(Time, Grd, id_pot_rho_3, wrk1(:,:,:)) endif endif if (id_pot_rho_4 > 0) then if (ind_potrho_4 .gt. 0) then call diagnose_3d(Time, Grd, id_pot_rho_4, T_diag(ind_potrho_4)%field(:,:,:)) else wrk1(:,:,:) = potential_density(Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), 4000.0) call diagnose_3d(Time, Grd, id_pot_rho_4, wrk1(:,:,:)) endif endif end subroutine ocean_density_diag ! NAME="ocean_density_diag" !####################################################################### ! ! ! ! ! If TEOS-10 is being used then we need to multiply Preformed Salinity and the ! Salinity factor to obtain absolute salinity for use in the TEOS10 EOS. ! ! If not using TEOS10 EOS, then copy the practical salinity into the ! density_salinity field for use in the preTEOS10 EOS or the linear EOS. ! ! Note that halo values are not generally valid for taup1 until the halos ! for index_salt and index_Fdelta have been updated inside of ocean_model.F90. ! ! subroutine update_ocean_density_salinity(T_prog, time_level, Dens) type(ocean_prog_tracer_type), intent(in) :: T_prog(:) integer, intent(in) :: time_level type(ocean_density_type), intent(inout) :: Dens integer i,j,k if (Dens%use_teos10) then do k=1,nk do j=jsd,jed do i=isd,ied Dens%rho_salinity(i,j,k,time_level) = T_prog(index_salt)%field(i,j,k,time_level)* & (1.0 + T_prog(index_Fdelta)%field(i,j,k,time_level)) enddo enddo enddo else do k=1,nk do j=jsd,jed do i=isd,ied Dens%rho_salinity(i,j,k,time_level) = T_prog(index_salt)%field(i,j,k,time_level) enddo enddo enddo endif end subroutine update_ocean_density_salinity ! NAME="update_ocean_density_salinity" !####################################################################### ! ! ! ! Diagnose pressure_at_depth and ocean density. ! Also send some diagnostics to diagnostic manager. ! ! subroutine update_ocean_density(Time, Thickness, Temp, Salt, Ext_mode, Dens, L_system, use_blobs) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Temp type(ocean_prog_tracer_type), intent(in) :: Salt type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_density_type), intent(inout) :: Dens type(ocean_lagrangian_type), intent(in) :: L_system logical, intent(in) :: use_blobs type(time_type) :: next_time integer :: i, j, k integer :: tau, taup1 integer :: stdoutunit stdoutunit=stdout() tau = Time%tau taup1 = Time%taup1 ! compute pressure at depth (dbars). ! do so consistently with all fields at time tau. if(vert_coordinate==GEOPOTENTIAL) then if (use_blobs) then Dens%pressure_at_depth(:,:,:) = pressure_in_dbars_blob(Thickness, Dens%rho(:,:,:,tau), & L_system%rho_dztup(:,:,:), L_system%rho_dztlo(:,:,:), & Ext_mode%patm_t(:,:,tau)+grav*Dens%rho(:,:,1,tau)*Ext_mode%eta_t(:,:,tau)) else Dens%pressure_at_depth(:,:,:) = pressure_in_dbars(Thickness, Dens%rho(:,:,:,tau), & Ext_mode%patm_t(:,:,tau)+grav*Dens%rho(:,:,1,tau)*Ext_mode%eta_t(:,:,tau)) endif elseif(vert_coordinate==ZSTAR .or. vert_coordinate==ZSIGMA) then if (use_blobs) then Dens%pressure_at_depth(:,:,:) = pressure_in_dbars_blob(Thickness, Dens%rho(:,:,:,tau), & L_system%rho_dztup(:,:,:), L_system%rho_dztlo(:,:,:), & Ext_mode%patm_t(:,:,tau)) else Dens%pressure_at_depth(:,:,:) = pressure_in_dbars(Thickness, Dens%rho(:,:,:,tau), & Ext_mode%patm_t(:,:,tau)) endif elseif(vert_coordinate==PRESSURE) then do k=1,nk do j=jsd,jed do i=isd,ied Dens%pressure_at_depth(i,j,k) = c2dbars*Thickness%depth_st(i,j,k) enddo enddo enddo elseif(vert_coordinate==PSTAR) then do k=1,nk do j=jsd,jed do i=isd,ied Dens%pressure_at_depth(i,j,k) = c2dbars*( & Ext_mode%patm_t(i,j,tau) & + Thickness%depth_st(i,j,k)*Thickness%pbot0r(i,j) & *(Ext_mode%pbot_t(i,j,tau)-Ext_mode%patm_t(i,j,tau)) & ) enddo enddo enddo elseif(vert_coordinate==PSIGMA) then do k=1,nk do j=jsd,jed do i=isd,ied Dens%pressure_at_depth(i,j,k) = c2dbars*( & Ext_mode%patm_t(i,j,tau) & + Thickness%depth_st(i,j,k) & *(Ext_mode%pbot_t(i,j,tau)-Ext_mode%patm_t(i,j,tau)) & ) enddo enddo enddo endif ! diagnose in situ density Dens%rho(:,:,:,taup1) = Grd%tmask(:,:,:) & *density(Dens%rho_salinity(:,:,:,taup1), Temp%field(:,:,:,taup1), Dens%pressure_at_depth(:,:,:)) if(rho0_density) then Dens%rho(:,:,:,taup1) = Grd%tmask(:,:,:)*rho0 endif ! compute drhodT and drhodS and drhodP ! do so with tau values for salinity, temperature, and pressure call density_derivs(Dens%rho(:,:,:,tau), Dens%rho_salinity(:,:,:,tau), & Temp%field(:,:,:,tau), Dens%pressure_at_depth(:,:,:), & Time, Dens%drhodT(:,:,:), Dens%drhodS(:,:,:), Dens%drhodP(:,:,:)) ! compute dpotrhodT and dpotrhodS ! do so with tau values for salinity, temperature, and pressure call density_derivs_potrho(Time, Dens%potrho(:,:,:), Dens%rho_salinity(:,:,:,tau),& Temp%field(:,:,:,tau), Dens%dpotrhodT(:,:,:), & Dens%dpotrhodS(:,:,:), Dens%dpotrhodP(:,:,:)) ! compute buoyancy frequency for diagnostic purposes if (use_blobs) then call compute_buoyfreq(Time, Thickness, Salt%fieldT(:,:,:), Temp%fieldT(:,:,:), Dens, use_blobs) else ! compute buoyancy frequency for diagnostic purposes call compute_buoyfreq(Time, Thickness, Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), Dens, use_blobs) endif call compute_density_diagnostics(Time, Thickness, Temp, Ext_mode, Dens, use_blobs) call compute_drhodxy(Time, Dens%rho_salinity(:,:,:,tau), Temp%field(:,:,:,tau), Dens) end subroutine update_ocean_density ! NAME="update_ocean_density" !####################################################################### ! ! ! ! Compute density for all grid points. ! ! Note that pressure here is ! ! sea pressure = absolute pressure - press_standard (dbars) ! ! and salinity is in model units (psu or g/kg). ! ! ! function density_field (salinity, theta, press) real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(in) :: press real, dimension(isd:ied,jsd:jed,nk) :: density_field real, dimension(isd:ied,jsd:jed,nk) :: pressure integer :: i, j, k real :: t1, t2, s1, sp5, p1, p1t1 real :: num, den if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_field): module must be initialized') endif if(density_equal_potrho) then pressure(:,:,:) = Grd%tmask(:,:,:)*potrho_press else pressure(:,:,:) = press(:,:,:) endif if(eos_linear) then do k=1,nk do j=jsd,jed do i=isd,ied density_field(i,j,k) = rho0 - alpha_linear_eos*theta(i,j,k) + beta_linear_eos*salinity(i,j,k) enddo enddo enddo elseif(eos_preteos10) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 num = a0 + t1*(a1 + t1*(a2+a3*t1) ) & + s1*(a4 + a5*t1 + a6*s1) & + p1*(a7 + a8*t2 + a9*s1 + p1*(a10+a11*t2)) den = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) & + p1*(b10 + p1t1*(b11*t2 + b12*p1)) denominator_r(i,j,k) = 1.0/(epsln+den) density_field(i,j,k) = num*denominator_r(i,j,k) enddo enddo enddo else ! eos_teos10 do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 num = v01 + t1*(v02 + t1*(v03 + v04*t1)) & + s1*(v05 + t1*(v06 + v07*t1) & + sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) & + p1*(v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(v17 + t1*(v18 + v19*t1) + v20*s1)) den = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1))) & + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*(v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + p1*(v47 + v48*t1))) denominator_r(i,j,k) = 1.0/(epsln+den) density_field(i,j,k) = num*denominator_r(i,j,k) enddo enddo enddo endif end function density_field ! NAME="density_field" !####################################################################### ! ! ! ! Compute density at a particular k-level. ! ! Note that pressure here is ! ! sea pressure = absolute pressure - press_standard (dbars) ! ! ! function density_level(salinity, theta, press) real, dimension(isd:,jsd:), intent(in) :: salinity real, dimension(isd:,jsd:), intent(in) :: theta real, dimension(isd:,jsd:), intent(in) :: press real, dimension(isd:ied,jsd:jed) :: density_level real, dimension(isd:ied,jsd:jed) :: pressure integer :: i, j real :: t1, t2, s1, sp5, p1, p1t1 real :: num, den if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_level): module must be initialized') endif if(density_equal_potrho) then pressure(:,:) = Grd%tmask(:,:,1)*potrho_press else pressure(:,:) = press(:,:) endif if(eos_linear) then do j=jsd,jed do i=isd,ied density_level(i,j) = rho0 - alpha_linear_eos*theta(i,j) + beta_linear_eos*salinity(i,j) enddo enddo elseif(eos_preteos10) then do j=jsd,jed do i=isd,ied t1 = theta(i,j) t2 = t1*t1 s1 = salinity(i,j) sp5 = sqrt(s1) p1 = pressure(i,j) - press_standard p1t1 = p1*t1 num = a0 + t1*(a1 + t1*(a2+a3*t1) ) & + s1*(a4 + a5*t1 + a6*s1) & + p1*(a7 + a8*t2 + a9*s1 + p1*(a10+a11*t2)) den = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) & + p1*(b10 + p1t1*(b11*t2 + b12*p1)) density_level(i,j) = num/(epsln+den) enddo enddo else ! eos_teos10 do j=jsd,jed do i=isd,ied t1 = theta(i,j) t2 = t1*t1 s1 = salinity(i,j) sp5 = sqrt(s1) p1 = pressure(i,j) - press_standard p1t1 = p1*t1 num = v01 + t1*(v02 + t1*(v03 + v04*t1)) & + s1*(v05 + t1*(v06 + v07*t1) & + sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) & + p1*(v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(v17 + t1*(v18 + v19*t1) + v20*s1)) den = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1))) & + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*(v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + p1*(v47 + v48*t1))) density_level(i,j) = num/(epsln+den) enddo enddo endif end function density_level ! NAME="density_level" !####################################################################### ! ! ! ! Compute density at a particular k-level and j index. This scheme ! is used in the vectorized version of the full convection scheme. ! ! Note that pressure here is ! ! sea pressure = absolute pressure - press_standard ! ! ! function density_line(salinity, theta, press) real, dimension(isd:), intent(in) :: salinity real, dimension(isd:), intent(in) :: theta real, dimension(isd:), intent(in) :: press real, dimension(isd:ied) :: density_line real, dimension(isd:ied) :: pressure integer :: i real :: t1, t2, s1, sp5, p1, p1t1 real :: num, den if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_line): module must be initialized') endif if(density_equal_potrho) then pressure(:) = potrho_press else pressure(:) = press(:) endif if(eos_linear) then do i=isd,ied density_line(i) = rho0 - alpha_linear_eos*theta(i) + beta_linear_eos*salinity(i) enddo elseif(eos_preteos10) then do i=isd,ied t1 = theta(i) t2 = t1*t1 s1 = salinity(i) sp5 = sqrt(s1) p1 = pressure(i) - press_standard p1t1 = p1*t1 num = a0 + t1*(a1 + t1*(a2+a3*t1) ) & + s1*(a4 + a5*t1 + a6*s1) & + p1*(a7 + a8*t2 + a9*s1 + p1*(a10+a11*t2)) den = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) & + p1*(b10 + p1t1*(b11*t2 + b12*p1)) density_line(i) = num/(epsln+den) enddo else ! eos_teos10 do i=isd,ied t1 = theta(i) t2 = t1*t1 s1 = salinity(i) sp5 = sqrt(s1) p1 = pressure(i) - press_standard p1t1 = p1*t1 num = v01 + t1*(v02 + t1*(v03 + v04*t1)) & + s1*(v05 + t1*(v06 + v07*t1) & + sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) & + p1*(v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(v17 + t1*(v18 + v19*t1) + v20*s1)) den = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1))) & + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*(v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + p1*(v47 + v48*t1))) density_line(i) = num/(epsln+den) enddo endif end function density_line ! NAME="density_line" !####################################################################### ! ! ! ! Compute neutral density for use in various layer diagnostics. ! ! Two options are presently available: ! ! A/ use rational polynomial (to be done) ! ! B/ use potential density referenced to pressure potrho_press ! McDougall recommends potential density referenced to 2000dbar, ! since the rational polynomial is not too good. ! ! Note that presently, the rational polynomial method defaults to ! potential density referenced to 2000dbar. ! The polynomial approximation from McDougall and Jackett (2005) ! is not recommended (as per Trevor McDougall, 2011). A new polynomial ! is being constructed and should be ready end of 2011. ! ! ! function neutral_density_field (salinity, theta) real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:ied,jsd:jed,nk) :: neutral_density_field integer :: i, j, k if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (neutral_density_field): module must be initialized') endif if(neutral_density_theta) then do k=1,nk do j=jsd,jed do i=isd,ied neutral_density_field(i,j,k) = theta(i,j,k) enddo enddo enddo elseif(eos_linear) then do k=1,nk do j=jsd,jed do i=isd,ied neutral_density_field(i,j,k) = rho0 - alpha_linear_eos*theta(i,j,k) + beta_linear_eos*salinity(i,j,k) enddo enddo enddo else if (diff_nrho_press) then neutral_density_field = potential_density(salinity,theta,nrho_press) else neutral_density_field = potential_density(salinity,theta,potrho_press) endif endif end function neutral_density_field ! NAME="neutral_density_field" !####################################################################### ! ! ! ! Compute neutral density for use in various layer diagnostics. ! ! Only test here the rational polynomial ! approximation given by McDougall and Jackett (2005). ! This test needs to be updated. ! ! ! function neutral_density_point (salinity, theta) real, intent(in) :: salinity real, intent(in) :: theta real :: neutral_density_point real :: t1, t2, s1, sp5 real :: num, den if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (neutral_density_point): module must be initialized') endif if(neutral_density_theta) then neutral_density_point = theta elseif(eos_linear) then neutral_density_point = rho0 - alpha_linear_eos*theta + beta_linear_eos*salinity else t1 = theta t2 = t1*t1 s1 = salinity sp5 = sqrt(s1) num = a0n + t1*(a1n + t1*(a2n+a3n*t1) ) & + s1*(a4n + a5n*t1 + a6n*s1) den = b0n + t1*(b1n + t1*(b2n + t1*(b3n + t1*b4n))) & + s1*(b5n + t1*(b6n + b7n*t2) + sp5*(b8n + b9n*t2)) neutral_density_point = num/(epsln+den) endif end function neutral_density_point ! NAME="neutral_density_point" !####################################################################### ! ! ! ! Compute potential density referenced to some given sea pressure. ! ! Note that potential density referenced to the surface (i.e., sigma_0) ! has a zero sea pressure, so pressure=0.0 should be the argument ! to the function. ! ! Note that pressure here is ! sea pressure = absolute pressure - press_standard (dbars) ! ! input pressure < 0 is an error, and model is brought down. ! ! ! function potential_density (salinity, theta, pressure) real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, intent(in) :: pressure real, dimension(isd:ied,jsd:jed,nk) :: potential_density integer :: i, j, k real :: t1, t2, s1, sp5, p1, p1t1 real :: num, den if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (potential_density): module must be initialized') endif if(pressure < 0.0) then call mpp_error(FATAL, & '==>Error in ocean_density_mod: potential density at pressure < 0 is not defined') endif if(eos_linear) then do k=1,nk do j=jsd,jed do i=isd,ied potential_density(i,j,k) = rho0 - alpha_linear_eos*theta(i,j,k) + beta_linear_eos*salinity(i,j,k) enddo enddo enddo elseif(eos_preteos10) then p1 = pressure if(p1 > 0.0) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1t1 = p1*t1 num = a0 + t1*(a1 + t1*(a2+a3*t1) ) & + s1*(a4 + a5*t1 + a6*s1) & + p1*(a7 + a8*t2 + a9*s1 + p1*(a10+a11*t2)) den = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) & + p1*(b10 + p1t1*(b11*t2 + b12*p1)) potential_density(i,j,k) = num/(epsln+den) enddo enddo enddo elseif(p1==0.0) then ! simplification for sigma_0 do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) num = a0 + t1*(a1 + t1*(a2+a3*t1) ) & + s1*(a4 + a5*t1 + a6*s1) den = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) potential_density(i,j,k) = num/(epsln+den) enddo enddo enddo endif ! endif for value of pressure else ! eos_teos10 p1 = pressure if(p1 > 0.0) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1t1 = p1*t1 num = v01 + t1*(v02 + t1*(v03 + v04*t1)) & + s1*(v05 + t1*(v06 + v07*t1) & + sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) & + p1*(v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(v17 + t1*(v18 + v19*t1) + v20*s1)) den = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1))) & + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*(v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + p1*(v47 + v48*t1))) potential_density(i,j,k) = num/(epsln+den) enddo enddo enddo elseif(p1==0.0) then ! for sigma_0 do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) num = v01 + t1*(v02 + t1*(v03 + v04*t1)) & + s1*(v05 + t1*(v06 + v07*t1) & + sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) den = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1))) & + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) potential_density(i,j,k) = num/(epsln+den) enddo enddo enddo endif ! endif for pressure endif ! endif for eos choice end function potential_density ! NAME="potential_density" !####################################################################### ! ! ! ! Diagnostics related to density. ! ! subroutine compute_density_diagnostics(Time, Thickness, Temp, Ext_mode, Dens, use_blobs) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Temp type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_density_type), intent(in) :: Dens logical, intent(in) :: use_blobs type(time_type) :: next_time logical :: used integer :: i, j, k integer :: tau, taup1 real :: mass, massip1, massjp1 real :: density_tau =0.0 real :: density_taup1 =0.0 real :: volume_tau =0.0 real :: volume_taup1 =0.0 real :: volume_taup12 =0.0 real :: mass_tau =0.0 real :: mass_taup1 =0.0 real :: pbot_adjust =0.0 real :: eta_adjust =0.0 real :: eta_adjust_approx =0.0 real :: eta_adjust_cstvolume=0.0 integer :: stdoutunit stdoutunit=stdout() tau = Time%tau taup1 = Time%taup1 ! compute cabbeling and thermobaricity parameters for diagnostics if(id_cabbeling > 0 .or. id_thermobaricity > 0) then call cabbeling_thermobaricity(Dens%rho(:,:,:,tau), Dens%rho_salinity(:,:,:,tau), & Temp%field(:,:,:,tau), Dens%pressure_at_depth(:,:,:), & Time, Dens%drhodT(:,:,:), Dens%drhodS(:,:,:), Dens%drhodP(:,:,:)) endif call diagnose_3d(Time, Grd, id_press, Dens%pressure_at_depth(:,:,:)) if (use_blobs) then call diagnose_3d(Time, Grd, id_rho, Dens%rhoT(:,:,:)) call diagnose_3d(Time, Grd, id_rhoE, Dens%rho(:,:,:,tau)) else call diagnose_3d(Time, Grd, id_rho, Dens%rho(:,:,:,tau)) endif call diagnose_3d(Time, Grd, id_eos_salinity, Dens%rho_salinity(:,:,:,tau)) if(id_int_rhodz > 0) then wrk1_2d(:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec wrk1_2d(i,j) = wrk1_2d(i,j) + Thickness%dzt(i,j,k)*Dens%rho(i,j,k,tau) enddo enddo enddo call diagnose_2d(Time, Grd, id_int_rhodz, wrk1_2d(:,:)) endif next_time = increment_time(Time%model_time, int(dtts), 0) if (need_data(id_pot_rho_wt,next_time)) then wrk1(:,:,:) = 0.0 do j=jsc,jec do i=isc,iec wrk1(i,j,nk) = Dens%potrho(i,j,nk) enddo enddo do k=1,nk-1 do j=jsc,jec do i=isc,iec if(Grd%tmask(i,j,k) > 0.0) then wrk1(i,j,k) = ( Grd%tmask(i,j,k+1)*Dens%potrho(i,j,k+1)*Thickness%rho_dzt(i,j,k+1,tau) & +Dens%potrho(i,j,k) *Thickness%rho_dzt(i,j,k,tau)) & /( epsln + Thickness%rho_dzt(i,j,k,tau) + Grd%tmask(i,j,k+1)*Thickness%rho_dzt(i,j,k+1,tau) ) endif enddo enddo enddo call diagnose_3d(Time, Grd, id_pot_rho_wt, wrk1(:,:,:)) endif if (need_data(id_pot_rho_et,next_time)) then wrk1(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec mass = Thickness%rho_dzt(i,j,k,tau) *Grd%dat(i,j) *Grd%tmask(i,j,k) massip1 = Thickness%rho_dzt(i+1,j,k,tau)*Grd%dat(i+1,j)*Grd%tmask(i+1,j,k) wrk1(i,j,k) = ( Dens%potrho(i,j,k)*mass + Dens%potrho(i+1,j,k)*massip1 ) & /( epsln + mass + massip1 ) enddo enddo enddo call diagnose_3d(Time, Grd, id_pot_rho_et, wrk1(:,:,:)) endif if (need_data(id_pot_rho_nt,next_time)) then wrk1(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec mass = Thickness%rho_dzt(i,j,k,tau) *Grd%dat(i,j) *Grd%tmask(i,j,k) massjp1 = Thickness%rho_dzt(i,j+1,k,tau)*Grd%dat(i,j+1)*Grd%tmask(i,j+1,k) wrk1(i,j,k) = ( Dens%potrho(i,j,k)*mass + Dens%potrho(i,j+1,k)*massjp1 ) & /( epsln + mass + massjp1 ) enddo enddo enddo call diagnose_3d(Time, Grd, id_pot_rho_nt, wrk1(:,:,:)) endif if (id_pbot_adjust>0 .or. id_eta_adjust>0 .or. & id_eta_adjust_approx>0 .or. id_eta_adjust_cstvolume>0 .or. id_rhoave>0) then ! These diagnostics are used to adjust the bottom pressure and surface height computed ! from a Boussinesq model. The adjustments correct, globally, for the spurious ! mass sources appearing in the Boussinesq fluid and the missing steric effect. ! The adjustments are NOT needed with non-Boussinesq models (pressure based vert coord), ! since the non-Boussinesq fluid correctly computes the bottom pressure and ! surface height according to mass conserving kinematics. ! ! In a Bouss model (depth based vert coordinates), when global mean density increases, ! so does the Boussinesq mass. So the pbot_adjust is negative in this case, in ! order to globally counteract effects from the spurious mass source. ! ! In a non-Bouss model, when global mean density increases, global mean sea level decreases. ! This global steric effect is absent from the Boussinesq fluid. To incorporate this ! this missing contribution to sea level, eta_adjust will be negative in this case. ! ! note that dzt has already been updated to its taup1 value. that is why we use ! rho0r*rho_dzt = dzt (when depth-based vertical coordinates). mass_tau = 0.0 mass_taup1 = 0.0 volume_tau = 0.0 volume_taup1 = 0.0 volume_taup12 = 0.0 rhodz_tau(:,:) = 0.0 rhodz_taup1(:,:) = 0.0 wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 wrk3_2d(:,:) = 0.0 wrk4_2d(:,:) = 0.0 ! do the vertical sum of rho*dzt if(vert_coordinate==GEOPOTENTIAL .or. vert_coordinate==ZSTAR .or. vert_coordinate==ZSIGMA) then do k=1,nk do j=jsc,jec do i=isc,iec rhodz_tau(i,j) = rhodz_tau(i,j) + Grd%tmask(i,j,k)*rho0r*Thickness%rho_dzt(i,j,k,tau) & *Dens%rho(i,j,k,tau) rhodz_taup1(i,j)= rhodz_taup1(i,j)+ Grd%tmask(i,j,k)*rho0r*Thickness%rho_dzt(i,j,k,taup1) & *Dens%rho(i,j,k,taup1) enddo enddo enddo else k=1 do j=jsc,jec do i=isc,iec rhodz_tau(i,j) = Grd%tmask(i,j,k)*grav_r*(Ext_mode%pbot_t(i,j,tau) -Ext_mode%patm_t(i,j,tau)) rhodz_taup1(i,j) = Grd%tmask(i,j,k)*grav_r*(Ext_mode%pbot_t(i,j,taup1)-Ext_mode%patm_t(i,j,taup1)) enddo enddo endif ! compute the mass and volume of a seawater column do j=jsd,jed do i=isd,ied wrk1_2d(i,j) = Grd%tmask(i,j,1)*Grd%dat(i,j)*rhodz_tau(i,j) wrk2_2d(i,j) = Grd%tmask(i,j,1)*Grd%dat(i,j)*rhodz_taup1(i,j) wrk3_2d(i,j) = Grd%tmask(i,j,1)*Grd%dat(i,j)*(Grd%ht(i,j) + Ext_mode%eta_t(i,j,tau)) wrk4_2d(i,j) = Grd%tmask(i,j,1)*Grd%dat(i,j)*(Grd%ht(i,j) + Ext_mode%eta_t(i,j,taup1)) enddo enddo ! perform global sums to get global mass and volume for ocean mass_tau = mpp_global_sum(Dom%domain2d,wrk1_2d(:,:),NON_BITWISE_EXACT_SUM) volume_tau = mpp_global_sum(Dom%domain2d,wrk3_2d(:,:),NON_BITWISE_EXACT_SUM) density_tau = mass_tau/(epsln+volume_tau) mass_taup1 = mpp_global_sum(Dom%domain2d,wrk2_2d(:,:),NON_BITWISE_EXACT_SUM) volume_taup1 = mpp_global_sum(Dom%domain2d,wrk4_2d(:,:),NON_BITWISE_EXACT_SUM) density_taup1 = mass_taup1/(epsln+volume_taup1) volume_taup12 = 0.5*(volume_taup1+volume_tau) if(id_pbot_adjust > 0) then pbot_adjust = -grav*c2dbars*area_total_r*volume_taup1*(density_taup1-density_tau) used = send_data(id_pbot_adjust, pbot_adjust, Time%model_time) endif if(id_eta_adjust > 0) then eta_adjust = area_total_r*volume_taup12*log((density_tau+epsln)/(density_taup1+epsln)) used = send_data(id_eta_adjust, eta_adjust, Time%model_time) endif if(id_eta_adjust_approx > 0) then eta_adjust_approx = area_total_r*volume_taup1*(density_tau/(epsln+density_taup1)-1.0) used = send_data(id_eta_adjust_approx, eta_adjust_approx, Time%model_time) endif if(id_eta_adjust_cstvolume > 0) then eta_adjust_cstvolume = area_total_r*Grd%tcellv*(density_tau/(epsln+density_taup1)-1.0) used = send_data(id_eta_adjust_cstvolume, eta_adjust_cstvolume, Time%model_time) endif if(id_rhoave > 0) then used = send_data(id_rhoave, density_tau, Time%model_time) endif endif if(debug_this_module) then write(stdoutunit,'(a)') ' ' write(stdoutunit,*) 'From ocean_density_mod: density chksums from compute_density_diagnostics' call write_timestamp(Time%model_time) call ocean_density_chksum(Time, Dens, use_blobs) endif end subroutine compute_density_diagnostics ! NAME="compute_density_diagnostics" !####################################################################### ! ! ! ! ! 1/ Compute ratio |grad neutral rho| / |grad local ref pot rho| ! for use in tform water mass analysis as per Iudicone et al. (2008). ! ! 2/ Compute rho*Area(h)/gamma_h, where "h" is the direction where ! gamma has the largest stratification, and where gamma is the ! locally referenced potential density. ! ! ! subroutine compute_diagnostic_factors (Time, Thickness, Dens, salinity, theta) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(inout) :: Dens real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta integer :: i, j, k, kp1, km1, tau real :: tmp, active_cells if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (compute_diagnostic_factors): module must be initialized') endif tau = Time%tau if (.not. update_diagnostic_factors) then do k=1,nk do j=jsd,jed do i=isd,ied Dens%watermass_factor(i,j,k) = neutralrho_interval_r*Grd%tmask(i,j,k) Dens%stratification_factor(i,j,k) = Grd%dxt(i,j)*Grd%dyt(i,j)*Thickness%dzt(i,j,k) Dens%rho_dztr_tau(i,j,k) = Thickness%rho_dztr(i,j,k) enddo enddo enddo return endif ! based on thinnest cell in a column, which is typically at k=1 or partial bottom if(smax_min_in_column) then do k=1,nk do j=jsd,jed do i=isd,ied smax(i,j,k) = Grd%tmask(i,j,k)*min(Thickness%dzt(i,j,k),Thickness%dzt(i,j,1))*dhorz_r(i,j) enddo enddo enddo ! straightforward aspect ratio calculation (generally the recommended method) else do k=1,nk do j=jsd,jed do i=isd,ied smax(i,j,k) = Grd%tmask(i,j,k)*Thickness%dzt(i,j,k)*dhorz_r(i,j) enddo enddo enddo endif ! over-ride the above using a prescribed smax if(smax_diag > 0.0) then do k=1,nk do j=jsd,jed do i=isd,ied smax(i,j,k) = min(1.0,smax_diag) enddo enddo enddo endif call diagnose_3d(Time, Grd, id_smax_dianeutral, smax(:,:,:)) ! some initialization do k=1,nk do j=jsd,jed do i=isd,ied wrk1_v(i,j,k,1) = 0.0 wrk1_v(i,j,k,2) = 0.0 wrk2_v(i,j,k,1) = 0.0 wrk2_v(i,j,k,2) = 0.0 wrk3_v(i,j,k,1) = 0.0 wrk3_v(i,j,k,2) = 0.0 wrk4_v(i,j,k,1) = 0.0 wrk4_v(i,j,k,2) = 0.0 wrk1(i,j,k) = 0.0 wrk2(i,j,k) = 0.0 wrk3(i,j,k) = 0.0 Dens%watermass_factor(i,j,k) = neutralrho_interval_r*Grd%tmask(i,j,k) wrk4(i,j,k) = Grd%tmask(i,j,k)/Thickness%dzwt(i,j,k) Dens%rho_dztr_tau(i,j,k) = Thickness%rho_dztr(i,j,k) enddo enddo enddo ! spatial derivatives of locally referenced potential density do k=1,nk kp1 = min(k+1,nk) wrk1_v(:,:,k,1) = FDX_T(theta(:,:,k)) *FMX(Grd%tmask(:,:,k)) wrk1_v(:,:,k,2) = FDX_T(salinity(:,:,k))*FMX(Grd%tmask(:,:,k)) wrk2_v(:,:,k,1) = FDY_T(theta(:,:,k)) *FMY(Grd%tmask(:,:,k)) wrk2_v(:,:,k,2) = FDY_T(salinity(:,:,k))*FMY(Grd%tmask(:,:,k)) do j=jsc,jec do i=isc,iec wrk3_v(i,j,k,1) = Grd%tmask(i,j,kp1)*wrk4(i,j,k)*(theta(i,j,k)-theta(i,j,kp1)) wrk3_v(i,j,k,2) = Grd%tmask(i,j,kp1)*wrk4(i,j,k)*(salinity(i,j,k)-salinity(i,j,kp1)) wrk1(i,j,k) = abs( Dens%drhodT(i,j,k)*wrk1_v(i,j,k,1) + Dens%drhodS(i,j,k)*wrk1_v(i,j,k,2) ) wrk2(i,j,k) = abs( Dens%drhodT(i,j,k)*wrk2_v(i,j,k,1) + Dens%drhodS(i,j,k)*wrk2_v(i,j,k,2) ) wrk3(i,j,k) = abs( Dens%drhodT(i,j,k)*wrk3_v(i,j,k,1) + Dens%drhodS(i,j,k)*wrk3_v(i,j,k,2) ) wrk4_v(i,j,k,2) = sqrt(wrk1(i,j,k)**2 + wrk2(i,j,k)**2 + wrk3(i,j,k)**2) + epsln_drhodz enddo enddo enddo ! compute |grad neutral density| wrk5(:,:,:) = 0.0 wrk6(:,:,:) = 0.0 do k=1,nk kp1 = min(k+1,nk) wrk5(:,:,k) = FDX_T(Dens%neutralrho(:,:,k))*FMX(Grd%tmask(:,:,k)) wrk6(:,:,k) = FDY_T(Dens%neutralrho(:,:,k))*FMY(Grd%tmask(:,:,k)) do j=jsc,jec do i=isc,iec tmp = Grd%tmask(i,j,kp1)*wrk4(i,j,k)*(Dens%neutralrho(i,j,k)-Dens%neutralrho(i,j,kp1)) wrk4_v(i,j,k,1) = sqrt(wrk5(i,j,k)**2 + wrk6(i,j,k)**2 + tmp**2) enddo enddo enddo ! compute watermass factor if(.not. eos_linear .and. grad_nrho_lrpotrho_compute) then wrk5(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec wrk5(i,j,k) = max(grad_nrho_lrpotrho_min, & min(grad_nrho_lrpotrho_max,Grd%tmask(i,j,k)*wrk4_v(i,j,k,1)/wrk4_v(i,j,k,2))) Dens%watermass_factor(i,j,k) = wrk5(i,j,k)*Dens%watermass_factor(i,j,k) enddo enddo enddo endif call diagnose_3d(Time, Grd, id_grad_nrho, wrk4_v(:,:,:,1)) call diagnose_3d(Time, Grd, id_grad_lrpotrho, wrk4_v(:,:,:,2)) call diagnose_3d(Time, Grd, id_grad_nrho_lrpotrho, wrk5(:,:,:)) call diagnose_3d(Time, Grd, id_watermass_factor, Dens%watermass_factor(:,:,:)) ! For stratification_factor, move derivatives to centre ! of a grid cell through 2-point averaging. call mpp_update_domains (wrk1(:,:,:), Dom%domain2d, complete=.false.) call mpp_update_domains (wrk2(:,:,:), Dom%domain2d, complete=.true.) wrk4(:,:,:) = 0.0 wrk5(:,:,:) = 0.0 wrk6(:,:,:) = 0.0 do k=1,nk kp1 = min(k+1,nk) km1 = max(k-1,1) do j=jsc,jec do i=isc,iec active_cells = epsln + Grd%tmask(i-1,j,k) + Grd%tmask(i+1,j,k) wrk4(i,j,k) = Grd%tmask(i,j,k)*(wrk1(i-1,j,k)+wrk1(i,j,k))/active_cells active_cells = epsln + Grd%tmask(i,j-1,k) + Grd%tmask(i,j+1,k) wrk5(i,j,k) = Grd%tmask(i,j,k)*(wrk2(i,j-1,k)+wrk2(i,j,k))/active_cells active_cells = epsln + Grd%tmask(i,j,km1) + Grd%tmask(i,j,kp1) wrk6(i,j,k) = Grd%tmask(i,j,k)*(wrk3(i,j,km1)+wrk3(i,j,k))/active_cells enddo enddo enddo ! smoothing may be useful in some cases if(smooth_stratification_factor) then call mpp_update_domains (wrk4(:,:,:), Dom%domain2d, complete=.false.) call mpp_update_domains (wrk5(:,:,:), Dom%domain2d, complete=.false.) call mpp_update_domains (wrk6(:,:,:), Dom%domain2d, complete=.true.) do k=1,nk wrk4(:,:,k) = S2D(wrk4(:,:,k)) wrk5(:,:,k) = S2D(wrk5(:,:,k)) wrk6(:,:,k) = S2D(wrk6(:,:,k)) enddo endif ! compute stratification factor. ! Vertical stratification is scaled down by smax in the if-tests ! to account for grid aspect ratio, and/or specific setting for smax_diag. wrk1(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec if(wrk4(i,j,k) >= wrk5(i,j,k) .and. wrk4(i,j,k) >= smax(i,j,k)*wrk6(i,j,k) .and. wrk4(i,j,k) > 0.0) then Dens%stratification_factor(i,j,k) = Dens%rho(i,j,k,tau)*Grd%dyt(i,j)*Thickness%dzt(i,j,k)/wrk4(i,j,k) wrk1(i,j,k) = 1.0 elseif(wrk5(i,j,k) >= smax(i,j,k)*wrk6(i,j,k) .and. wrk5(i,j,k) >= wrk4(i,j,k) .and. wrk5(i,j,k) > 0.0) then Dens%stratification_factor(i,j,k) = Dens%rho(i,j,k,tau)*Grd%dxt(i,j)*Thickness%dzt(i,j,k)/wrk5(i,j,k) wrk1(i,j,k) = 2.0 elseif(smax(i,j,k)*wrk6(i,j,k) >= wrk4(i,j,k) .and. smax(i,j,k)*wrk6(i,j,k) >= wrk5(i,j,k) .and. wrk6(i,j,k) > 0.0) then Dens%stratification_factor(i,j,k) = Dens%rho(i,j,k,tau)*Grd%dxt(i,j)*Grd%dyt(i,j)/wrk6(i,j,k) wrk1(i,j,k) = 3.0 else Dens%stratification_factor(i,j,k) = 0.0 wrk1(i,j,k) = 0.0 endif enddo enddo enddo call diagnose_3d(Time, Grd, id_stratification_factor, Dens%stratification_factor(:,:,:)) call diagnose_3d(TIme, Grd, id_stratification_axis, wrk1(:,:,:)) end subroutine compute_diagnostic_factors ! NAME="compute_diagnostic_factors" !####################################################################### ! ! ! ! Compute density as a function of surface salinity, surface theta, ! and in situ sea pressure. ! ! Note that pressure here is ! sea pressure = absolute pressure - press_standard (dbars) ! ! For use in KPP mixed layer scheme ! ! function density_sfc (salinity, theta, press) real, intent(in), dimension(isd:,jsd:,:) :: salinity real, intent(in), dimension(isd:,jsd:,:) :: theta real, intent(in), dimension(isd:,jsd:,:) :: press real, dimension(isd:ied,jsd:jed,nk) :: density_sfc real, dimension(isd:ied,jsd:jed,nk) :: pressure integer :: i, j, k real :: t1, t2, s1, sp5, p1, p1t1 real :: num, den if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_sfc): module must be initialized') endif if(density_equal_potrho) then pressure(:,:,:) = Grd%tmask(:,:,:)*potrho_press else pressure(:,:,:) = press(:,:,:) endif if(eos_linear) then do k=1,nk do j=jsd,jed do i=isd,ied density_sfc(i,j,k) = rho0 - alpha_linear_eos*theta(i,j,1) + beta_linear_eos*salinity(i,j,1) enddo enddo enddo elseif(eos_preteos10) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,1) t2 = t1*t1 s1 = salinity(i,j,1) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 num = a0 + t1*(a1 + t1*(a2+a3*t1) ) & + s1*(a4 + a5*t1 + a6*s1) & + p1*(a7 + a8*t2 + a9*s1 + p1*(a10+a11*t2)) den = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) & + p1*(b10 + p1t1*(b11*t2 + b12*p1)) density_sfc(i,j,k) = num/(epsln + den) enddo enddo enddo else ! eos_teos10 do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,1) t2 = t1*t1 s1 = salinity(i,j,1) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 num = v01 + t1*(v02 + t1*(v03 + v04*t1)) & + s1*(v05 + t1*(v06 + v07*t1) & + sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) & + p1*(v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(v17 + t1*(v18 + v19*t1) + v20*s1)) den = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1))) & + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*(v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + p1*(v47 + v48*t1))) density_sfc(i,j,k) = num/(epsln + den) enddo enddo enddo endif end function density_sfc ! NAME="density_sfc" !####################################################################### ! ! ! ! Compute density at a single model grid point. ! ! Note that pressure here is ! ! sea pressure = absolute pressure - press_standard (dbars) ! ! ! function density_point (s1, t1, p1_dbars) real, intent(in) :: s1 real, intent(in) :: t1 real, intent(in) :: p1_dbars real :: t2, sp5, p1, p1t1 real :: num, den real :: density_point real :: pressure if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_point): module must be initialized') endif if(density_equal_potrho) then pressure = potrho_press else pressure = p1_dbars endif if(eos_linear) then density_point = rho0 - alpha_linear_eos*t1 + beta_linear_eos*s1 elseif(eos_preteos10) then t2 = t1*t1 sp5 = sqrt(s1) p1 = pressure - press_standard p1t1 = p1*t1 num = a0 + t1*(a1 + t1*(a2+a3*t1)) & + s1*(a4 + a5*t1 + a6*s1) & + p1*(a7 + a8*t2 + a9*s1 + p1*(a10+a11*t2)) den = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4 ))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) & + p1*(b10 + p1t1*(b11*t2 + b12*p1)) density_point = num/(epsln+den) else ! eos_teos10 t2 = t1*t1 sp5 = sqrt(s1) p1 = pressure - press_standard p1t1 = p1*t1 num = v01 + t1*(v02 + t1*(v03 + v04*t1)) & + s1*(v05 + t1*(v06 + v07*t1) & + sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) & + p1*(v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(v17 + t1*(v18 + v19*t1) + v20*s1)) den = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1))) & + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*(v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + p1*(v47 + v48*t1))) density_point = num/(epsln+den) endif end function density_point ! NAME="density_point" !####################################################################### ! ! ! ! Compute partial derivative of density with respect to ! temperature and with respect to salinity. Hold pressure constant. ! ! Pressure here is ! ! sea pressure = absolute press - press_standard (dbars) ! ! ! subroutine density_derivs_field (rho, salinity, theta, press, Time, & density_theta, density_salinity, density_press) type(ocean_time_type), intent(in) :: Time real, dimension(isd:,jsd:,:), intent(in) :: rho real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(in) :: press real, dimension(isd:,jsd:,:), intent(inout) :: density_theta real, dimension(isd:,jsd:,:), intent(inout) :: density_salinity real, dimension(isd:,jsd:,:), intent(inout) :: density_press real, dimension(isd:ied,jsd:jed,nk) :: pressure integer :: i, j, k real :: t1, t2, s1, sp5, p1, p1t1 real :: dnum_dtheta, dden_dtheta real :: dnum_dsalinity, dden_dsalinity real :: dnum_dpress, dden_dpress if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_derivs_field): module must be initialized') endif if(density_equal_potrho) then pressure(:,:,:) = Grd%tmask(:,:,:)*potrho_press else pressure(:,:,:) = press(:,:,:) endif if(eos_linear) then do k=1,nk do j=jsd,jed do i=isd,ied density_theta(i,j,k) = -alpha_linear_eos density_salinity(i,j,k) = beta_linear_eos density_press(i,j,k) = 0.0 enddo enddo enddo elseif(eos_preteos10) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 dnum_dtheta = a1 + t1*(two_a2 + three_a3*t1) & + a5*s1 & + p1t1*(two_a8 + two_a11*p1) dden_dtheta = b1 + t1*(two_b2 + t1*(three_b3 + four_b4*t1)) & + s1*(b6 + t1*(three_b7*t1 + two_b9*sp5)) & + p1*p1*(three_b11*t2 + b12*p1) dnum_dsalinity = a4 + a5*t1 + two_a6*s1 + a9*p1 dden_dsalinity = b5 + t1*(b6 + b7*t2) + sp5*(onep5_b8 + onep5_b9*t2) dnum_dpress = a7 + a9*s1 + two_a10*p1 + t2*(a8 + two_a11*p1) dden_dpress = b10 + p1*t1*(two_b11*t2 + three_b12*p1) density_theta(i,j,k) = denominator_r(i,j,k)*(dnum_dtheta - rho(i,j,k)*dden_dtheta) density_salinity(i,j,k) = denominator_r(i,j,k)*(dnum_dsalinity - rho(i,j,k)*dden_dsalinity) density_press(i,j,k) = denominator_r(i,j,k)*(dnum_dpress - rho(i,j,k)*dden_dpress) enddo enddo enddo else ! eos_teos10 do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 dnum_dtheta = v02 + t1*(two_v03 + three_v04*t1) & + s1*((v06+two_v07*t1) & + sp5*(v09 + t1*(two_v10 + three_v11*t1))) & + p1*((v13 + two_v14*t1) + s1*v16 & + p1*(v18+two_v19*t1)) dden_dtheta = v22 + t1*(two_v23 +t1*(three_v24+four_v25*t1)) & + s1*(v27 + t1*(two_v28 + t1*(three_v29 + four_v30*t1)) & + sp5*(v32+t1*(two_v33 + t1*(three_v34 + four_v35*t1)))) & + p1*((v38 + t1*(two_v39 + three_v40*t1)) & + s1*v42 & + p1*((v44 + two_v45*t1 + v46*s1) & + p1*v48)) dnum_dsalinity = & + (v05 + t1*(v06 + v07*t1) & + 1.5*sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1))))& + p1*((v15 + v16*t1) & + p1*v20) dden_dsalinity = & (v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + two_v36*s1 & + 1.5*sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*( v41 + v42*t1 + p1*t1*v46 ) dnum_dpress = v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(two_v17 + t1*(two_v18 + two_v19*t1) + two_v20*s1) dden_dpress = (v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(two_v43 + t1*(two_v44 + two_v45*t1 + two_v46*s1) & + p1*(three_v47 + three_v48*t1))) density_theta(i,j,k) = denominator_r(i,j,k)*(dnum_dtheta - rho(i,j,k)*dden_dtheta) density_salinity(i,j,k) = denominator_r(i,j,k)*(dnum_dsalinity - rho(i,j,k)*dden_dsalinity) density_press(i,j,k) = denominator_r(i,j,k)*(dnum_dpress - rho(i,j,k)*dden_dpress)*c2dbars enddo enddo enddo endif call diagnose_3d(TIme, Grd, id_drhodtheta, density_theta(:,:,:)) call diagnose_3d(Time, Grd, id_drhodsalt, density_salinity(:,:,:)) call diagnose_3d(Time, Grd, id_drhodpress, density_press(:,:,:)) if (id_sound_speed2 > 0) then wrk1(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec if(density_press(i,j,k) /= 0.0) then wrk1(i,j,k) = 1.0/(c2dbars*density_press(i,j,k)) endif enddo enddo enddo call diagnose_3d(Time, Grd, id_sound_speed2, wrk1(:,:,:)) endif if(id_thermal_expand > 0) then wrk1(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec wrk1(i,j,k) = -Grd%tmask(i,j,k)*density_theta(i,j,k)/(epsln+rho(i,j,k)) enddo enddo enddo call diagnose_3d(Time, Grd, id_thermal_expand, wrk1(:,:,:)) endif if(id_haline_contract > 0) then wrk1(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec wrk1(i,j,k) = Grd%tmask(i,j,k)*density_salinity(i,j,k)/(epsln+rho(i,j,k)) enddo enddo enddo call diagnose_3d(Time, Grd, id_haline_contract, wrk1(:,:,:)) endif end subroutine density_derivs_field ! NAME="density_derivs_field" !####################################################################### ! ! ! ! Compute partial derivative of density with respect to ! temperature and with respect to salinity. Use potrho_press ! reference depth for potential density rather than neutral. ! ! These coefficients are useful for certain diagnostics, such as APE ! calculations. ! ! The algorithm is identical to density_derivs_field, only here we set ! the pressure to a constant reference pressure value, potrho_press. ! ! Pressure here is ! sea pressure = absolute press - press_standard (dbars) ! ! ! subroutine density_derivs_potrho (Time, rho, salinity, theta, & density_theta, density_salinity, density_press) type(ocean_time_type), intent(in) :: Time real, dimension(isd:,jsd:,:), intent(in) :: rho real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(inout) :: density_theta real, dimension(isd:,jsd:,:), intent(inout) :: density_salinity real, dimension(isd:,jsd:,:), intent(inout) :: density_press real, dimension(isd:ied,jsd:jed,nk) :: pressure integer :: i, j, k real :: t1, t2, s1, sp5, p1, p1t1 real :: dnum_dtheta, dden_dtheta real :: dnum_dsalinity, dden_dsalinity real :: dnum_dpress, dden_dpress if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_derivs_field): module must be initialized') endif pressure(:,:,:) = Grd%tmask(:,:,:)*potrho_press if(eos_linear) then do k=1,nk do j=jsd,jed do i=isd,ied density_theta(i,j,k) = -alpha_linear_eos density_salinity(i,j,k) = beta_linear_eos density_press(i,j,k) = 0.0 enddo enddo enddo elseif(eos_preteos10) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 dnum_dtheta = a1 + t1*(two_a2 + three_a3*t1) & + a5*s1 & + p1t1*(two_a8 + two_a11*p1) dden_dtheta = b1 + t1*(two_b2 + t1*(three_b3 + four_b4*t1)) & + s1*(b6 + t1*(three_b7*t1 + two_b9*sp5)) & + p1*p1*(three_b11*t2 + b12*p1) dnum_dsalinity = a4 + a5*t1 + two_a6*s1 + a9*p1 dden_dsalinity = b5 + t1*(b6 + b7*t2) + sp5*(onep5_b8 + onep5_b9*t2) dnum_dpress = a7 + a9*s1 + two_a10*p1 + t2*(a8 + two_a11*p1) dden_dpress = b10 + p1*t1*(two_b11*t2 + three_b12*p1) density_theta(i,j,k) = denominator_r(i,j,k)*(dnum_dtheta - rho(i,j,k)*dden_dtheta) density_salinity(i,j,k) = denominator_r(i,j,k)*(dnum_dsalinity - rho(i,j,k)*dden_dsalinity) density_press(i,j,k) = denominator_r(i,j,k)*(dnum_dpress - rho(i,j,k)*dden_dpress) enddo enddo enddo else ! eos_teos10 do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = pressure(i,j,k) - press_standard p1t1 = p1*t1 dnum_dtheta = v02 + t1*(two_v03 + three_v04*t1) & + s1*((v06+two_v07*t1) & + sp5*(v09 + t1*(two_v10 + three_v11*t1))) & + p1*((v13 + two_v14*t1) + s1*v16 & + p1*(v18+two_v19*t1)) dden_dtheta = v22 + t1*(two_v23 +t1*(three_v24+four_v25*t1)) & + s1*(v27 + t1*(two_v28 + t1*(three_v29 + four_v30*t1)) & + sp5*(v32+t1*(two_v33 + t1*(three_v34 +four_v35*t1)))) & + p1*((v38 + t1*(two_v39 + three_v40*t1)) & + s1*v42 & + p1*((v44 + two_v45*t1 + v46*s1) & + p1*v48)) dnum_dsalinity = & + (v05 + t1*(v06 + v07*t1) & + 1.5*sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1))))& + p1*((v15 + v16*t1) & + p1*v20) dden_dsalinity = & (v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + two_v36*s1 & + 1.5*sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*( v41 + v42*t1 + p1*t1*v46 ) dnum_dpress = v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + p1*(two_v17 + t1*(two_v18 + two_v19*t1) + two_v20*s1) dden_dpress = (v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(two_v43 + t1*(two_v44 + two_v45*t1 + two_v46*s1) & + p1*(three_v47 + three_v48*t1))) density_theta(i,j,k) = denominator_r(i,j,k)*(dnum_dtheta - rho(i,j,k)*dden_dtheta) density_salinity(i,j,k) = denominator_r(i,j,k)*(dnum_dsalinity - rho(i,j,k)*dden_dsalinity) density_press(i,j,k) = denominator_r(i,j,k)*(dnum_dpress - rho(i,j,k)*dden_dpress)*c2dbars enddo enddo enddo endif call diagnose_3d(Time, Grd, id_dpotrhodtheta, density_theta(:,:,:)) call diagnose_3d(Time, Grd, id_dpotrhodsalt, density_salinity(:,:,:)) call diagnose_3d(Time, Grd, id_dpotrhodpress, density_press(:,:,:)) end subroutine density_derivs_potrho ! NAME="density_derivs_potrho" !####################################################################### ! ! ! ! Compute partial derivative of density with respect to ! temperature and with respect to salinity. Hold pressure constant. ! ! Pressure here is sea pressure = absolute press - press_standard ! ! ! subroutine density_derivs_level (rho, salinity, theta, press, Time, klevel, & density_theta, density_salinity, density_press) type(ocean_time_type), intent(in) :: Time real, dimension(isd:,jsd:), intent(in) :: rho real, dimension(isd:,jsd:), intent(in) :: salinity real, dimension(isd:,jsd:), intent(in) :: theta real, dimension(isd:,jsd:), intent(in) :: press integer, intent(in) :: klevel real, dimension(isd:,jsd:), intent(inout) :: density_theta real, dimension(isd:,jsd:), intent(inout) :: density_salinity real, dimension(isd:,jsd:), intent(inout) :: density_press real, dimension(isd:ied,jsd:jed) :: pressure integer :: i, j real :: t1, t2, s1, sp5, p1, p1t1 real :: dnum_dtheta, dden_dtheta real :: dnum_dsalinity, dden_dsalinity real :: dnum_dpress, dden_dpress if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_derivs_level): module must be initialized') endif if(density_equal_potrho) then pressure(:,:) = Grd%tmask(:,:,1)*potrho_press else pressure(:,:) = press(:,:) endif if(eos_linear) then do j=jsd,jed do i=isd,ied density_theta(i,j) = -alpha_linear_eos density_salinity(i,j) = beta_linear_eos density_press(i,j) = 0.0 enddo enddo elseif(eos_preteos10) then do j=jsd,jed do i=isd,ied t1 = theta(i,j) t2 = t1*t1 s1 = salinity(i,j) sp5 = sqrt(s1) p1 = pressure(i,j) - press_standard p1t1 = p1*t1 dnum_dtheta = a1 + t1*(two_a2 + three_a3*t1) & + a5*s1 & + p1t1*(two_a8 + two_a11*p1) dden_dtheta = b1 + t1*(two_b2 + t1*(three_b3 + four_b4*t1)) & + s1*(b6 + t1*(three_b7*t1 + two_b9*sp5)) & + p1*p1*(three_b11*t2 + b12*p1) dnum_dsalinity = a4 + a5*t1 + two_a6*s1 + a9*p1 dden_dsalinity = b5 + t1*(b6 + b7*t2) + sp5*(onep5_b8 + onep5_b9*t2) dnum_dpress = a7 + a9*s1 + two_a10*p1 + t2*(a8 + two_a11*p1) dden_dpress = b10 + p1*t1*(two_b11*t2 + three_b12*p1) density_theta(i,j) = denominator_r(i,j,klevel)*(dnum_dtheta - rho(i,j)*dden_dtheta) density_salinity(i,j) = denominator_r(i,j,klevel)*(dnum_dsalinity - rho(i,j)*dden_dsalinity) density_press(i,j) = denominator_r(i,j,klevel)*(dnum_dpress - rho(i,j)*dden_dpress) enddo enddo else ! eos_teos10 t1 = theta(i,j) t2 = t1*t1 s1 = salinity(i,j) sp5 = sqrt(s1) p1 = pressure(i,j) - press_standard p1t1 = p1*t1 dnum_dtheta = v02 + t1*(two_v03 + three_v04*t1) & + s1*((v06+two_v07*t1) & + sp5*(v09 + t1*(two_v10 + three_v11*t1))) & + p1*((v13 + two_v14*t1) + s1*v16 & + p1*(v18+two_v19*t1)) dden_dtheta = v22 + t1*(two_v23 +t1*(three_v24+four_v25*t1)) & + s1*(v27 + t1*(two_v28 + t1*(three_v29 + four_v30*t1)) & + sp5*(v32+t1*(two_v33 + t1*(three_v34 + four_v35*t1)))) & + p1*((v38 + t1*(two_v39 + three_v40*t1)) & + s1*v42 & + p1*((v44 + two_v45*t1 + v46*s1) & + p1*(v48))) dnum_dsalinity = & + (v05 + t1*(v06 + v07*t1) & + 1.5*sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1)))) & + p1*((v15 + v16*t1) & + p1*(v20)) dden_dsalinity = & (v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + two_v36*s1 & + 1.5*sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*( v41 + v42*t1 + p1*t1*v46 ) dnum_dpress = v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + 2.0*p1*(v17 + t1*(v18 + v19*t1) + v20*s1) dden_dpress = (v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + 2.*p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + 1.5*p1*(v47 + v48*t1))) density_theta(i,j) = denominator_r(i,j,klevel)*(dnum_dtheta - rho(i,j)*dden_dtheta) density_salinity(i,j) = denominator_r(i,j,klevel)*(dnum_dsalinity - rho(i,j)*dden_dsalinity) density_press(i,j) = denominator_r(i,j,klevel)*(dnum_dpress - rho(i,j)*dden_dpress) endif ! endif for eos choice end subroutine density_derivs_level ! NAME="density_derivs_level" !####################################################################### ! ! ! ! Compute partial derivative of density with respect to ! temperature and with respect to salinity. Do so here for a point. ! ! Pressure here is ! ! sea pressure = absolute pressure - press_standard (dbars) ! ! ! subroutine density_derivs_point (rho, salinity, theta, press,& density_theta, density_salinity, density_press) real, intent(in) :: rho real, intent(in) :: salinity real, intent(in) :: theta real, intent(in) :: press real, intent(out) :: density_theta real, intent(out) :: density_salinity real, intent(out) :: density_press real :: t1, t2, s1, sp5, p1, p1t1 real :: dnum_dtheta, dden_dtheta real :: dnum_dsalinity, dden_dsalinity real :: dnum_dpress, dden_dpress real :: denominator_point, denrecip real :: pressure if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_derivs_point): module must be initialized') endif if(density_equal_potrho) then pressure = potrho_press else pressure = press endif if(eos_linear) then density_theta = -alpha_linear_eos density_salinity = beta_linear_eos density_press = 0.0 elseif(eos_preteos10) then t1 = theta t2 = t1*t1 s1 = salinity sp5 = sqrt(s1) p1 = pressure - press_standard p1t1 = p1*t1 denominator_point = b0 + t1*(b1 + t1*(b2 + t1*(b3 + t1*b4 ))) & + s1*(b5 + t1*(b6 + b7*t2) + sp5*(b8 + b9*t2)) & + p1*(b10 + p1t1*(b11*t2 + b12*p1)) denrecip = 1.0/(denominator_point + epsln) dnum_dtheta = a1 + t1*(2.0*a2 + 3.0*a3*t1) & + a5*s1 & + 2.0*p1t1*(a8 + a11*p1) dden_dtheta = b1 + t1*(2.0*b2 + t1*(3.0*b3 + 4.0*b4*t1)) & + s1*(b6 + t1*(3.0*b7*t1 + 2.0*b9*sp5)) & + p1*p1*(3.0*b11*t2 + b12*p1) dnum_dsalinity = a4 + a5*t1 + two_a6*s1 + a9*p1 dden_dsalinity = b5 + t1*(b6 + b7*t2) + sp5*(onep5_b8 + onep5_b9*t2) dnum_dpress = a7 + a9*s1 + two_a10*p1 + t2*(a8 + two_a11*p1) dden_dpress = b10 + p1*t1*(two_b11*t2 + three_b12*p1) density_theta = denrecip*(dnum_dtheta - rho*dden_dtheta) density_salinity = denrecip*(dnum_dsalinity - rho*dden_dsalinity) density_press = denrecip*(dnum_dpress - rho*dden_dpress) else ! eos_teos10 t1 = theta t2 = t1*t1 s1 = salinity sp5 = sqrt(s1) p1 = pressure - press_standard p1t1 = p1*t1 denominator_point = v21 + t1*(v22 + t1*(v23 + t1*(v24 + v25*t1)))& + s1*(v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + v36*s1 & + sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*(v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + p1*(v47 + v48*t1))) denrecip = 1./(denominator_point + epsln ) dnum_dtheta = v02 + t1*(two_v03 + three_v04*t1) & + s1*(v06+two_v07*t1 & + sp5*(v09 + t1*(two_v10 + three_v11*t1))) & + p1*((v13 + two_v14*t1) + s1*v16 & + p1*(v18+two_v19*t1)) dden_dtheta = v22 + t1*(two_v23 +t1*(three_v24+four_v25*t1)) & + s1* (v27+t1*(two_v28 + t1*(three_v29 + four_v30*t1)) & + sp5*(v32+t1*(two_v33 + t1*(three_v34 + four_v35*t1)))) & + p1*(v38 + t1*(two_v39 + three_v40*t1)+ s1*v42 & + p1*(v44 + two_v45*t1 + v46*s1 & + p1*v48)) dnum_dsalinity = & + (v05 + t1*(v06 + v07*t1) & + 1.5*sp5*(v08 + t1*(v09 + t1*(v10 + v11*t1))))& + p1*((v15 + v16*t1) & + p1*(v20)) dden_dsalinity = & (v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + two_v36*s1 & + 1.5*sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*( v41 + v42*t1 + p1*t1*v46 ) dnum_dpress = v12 + t1*(v13 + v14*t1) + s1*(v15 + v16*t1) & + 2.0*p1*(v17 + t1*(v18 + v19*t1) + v20*s1) dden_dpress = (v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + 2.*p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + 1.5*p1*(v47 + v48*t1))) density_theta = denrecip*(dnum_dtheta - rho*dden_dtheta) density_salinity = denrecip*(dnum_dsalinity - rho*dden_dsalinity) density_press = denrecip*(dnum_dpress - rho*dden_dpress) endif ! endif for eos choice end subroutine density_derivs_point ! NAME="density_derivs_point" !####################################################################### ! ! ! ! Diagnostic sends for cabbeling and thermobaricity parameters. ! ! Pressure here is ! sea pressure = absolute press - press_standard (dbars) ! ! ! subroutine cabbeling_thermobaricity(rho, salinity, theta, press, Time, & density_theta, density_salinity, density_press) real, dimension(isd:,jsd:,:), intent(in) :: rho real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(in) :: press type(ocean_time_type), intent(in) :: Time real, dimension(isd:,jsd:,:), intent(in) :: density_theta real, dimension(isd:,jsd:,:), intent(in) :: density_salinity real, dimension(isd:,jsd:,:), intent(in) :: density_press real, dimension(isd:ied,jsd:jed,nk) :: cabbeling_param real, dimension(isd:ied,jsd:jed,nk) :: thermobaric_param real, dimension(isd:ied,jsd:jed,nk) :: pressure if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (cabbeling_thermobaricity): module must be initialized') endif if(density_equal_potrho) then pressure(:,:,:) = Grd%tmask(:,:,:)*potrho_press else pressure(:,:,:) = press(:,:,:) endif call calc_cabbeling_thermobaricity(rho, salinity, theta, pressure, & density_theta, density_salinity, density_press, cabbeling_param, thermobaric_param) call diagnose_3d(Time, Grd, id_cabbeling, cabbeling_param(:,:,:)) call diagnose_3d(Time, Grd, id_thermobaricity, thermobaric_param(:,:,:)) end subroutine cabbeling_thermobaricity ! NAME="cabbeling_thermobaricity" !####################################################################### ! ! ! ! Compute cabbeling and thermobaricity parameters, as defined in ! McDougall (1987). ! ! Pressure here is ! sea pressure = absolute press - press_standard (dbars) ! ! ! subroutine calc_cabbeling_thermobaricity(rho, salinity, theta, press, & density_theta, density_salinity, density_press, cabbeling_param, thermobaric_param) real, dimension(isd:,jsd:,:), intent(in) :: rho real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(in) :: press real, dimension(isd:,jsd:,:), intent(in) :: density_theta real, dimension(isd:,jsd:,:), intent(in) :: density_salinity real, dimension(isd:,jsd:,:), intent(in) :: density_press real, dimension(isd:,jsd:,:), intent(inout) :: cabbeling_param real, dimension(isd:,jsd:,:), intent(inout) :: thermobaric_param integer :: i, j, k real :: t1, t2, s1, sp5, spm5, p1, p2, p1t1 real :: rhotheta_rhosalinity real :: rho_inv real :: d2rho_dtheta2, d2rho_dsalin2, d2rho_dsalin_dtheta real :: d2rho_dsalin_dpress, d2rho_dtheta_dpress real :: dden_dtheta, dden_dsalinity, dden_dpress real :: d2num_dtheta2, d2num_dsalin2, d2num_dtheta_dsalin real :: d2num_dsalin_dpress, d2num_dtheta_dpress real :: d2den_dtheta2, d2den_dsalin2, d2den_dsalin_dtheta real :: d2den_dsalin_dpress, d2den_dtheta_dpress if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (calc_cabbeling_thermobaricity): module must be initialized') endif ! zero initialization, and values for when eos_linear is used cabbeling_param(:,:,:) = 0.0 ! cabbeling parameter thermobaric_param(:,:,:) = 0.0 ! thermobaricity parameter if(eos_preteos10) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = press(i,j,k) - press_standard p2 = p1*p1 p1t1 = p1*t1 if(s1 > 0.0) then spm5 = 1.0/sp5 else spm5 = 0.0 endif rhotheta_rhosalinity = Grd%tmask(i,j,k)*density_theta(i,j,k)/(epsln+density_salinity(i,j,k)) dden_dtheta = b1 + t1*(two_b2 + t1*(three_b3 + four_b4*t1)) & + s1*(b6 + t1*(three_b7*t1 + two_b9*sp5)) & + p1*p1*(three_b11*t2 + b12*p1) dden_dsalinity = b5 + t1*(b6 + b7*t2) + sp5*(onep5_b8 + onep5_b9*t2) dden_dpress = b10 + p1*t1*(two_b11*t2 + three_b12*p1) d2num_dtheta2 = two_a2 + six_a3*t1 + two_a8*p1 + two_a11*p2 d2num_dsalin2 = two_a6 d2num_dtheta_dsalin = a5 d2num_dsalin_dpress = a9 d2num_dtheta_dpress = two_a8*t1 + four_a11*p1*t1 d2den_dtheta2 = two_b2 + six_b3*t1 + twelve_b4*t2 & + six_b7*s1*t1 + two_b9*s1*sp5 + six_b11*p2*t1 d2den_dsalin2 = .75*spm5*(b8 + b9*t2) d2den_dsalin_dtheta = b6 + three_b7*t2 + three_b9*sp5*t1 d2den_dsalin_dpress = 0.0 d2den_dtheta_dpress = six_b11*p1*t2 + three_b12*p2 d2rho_dtheta2 = denominator_r(i,j,k) & *(d2num_dtheta2-2.0*density_theta(i,j,k)*dden_dtheta-rho(i,j,k)*d2den_dtheta2) d2rho_dsalin2 = denominator_r(i,j,k) & *(d2num_dsalin2-2.0*density_salinity(i,j,k)*dden_dsalinity-rho(i,j,k)*d2den_dsalin2) d2rho_dsalin_dtheta = denominator_r(i,j,k) & *( d2num_dtheta_dsalin & -density_salinity(i,j,k)*dden_dtheta & -density_theta(i,j,k)*dden_dsalinity & -rho(i,j,k)*d2den_dsalin_dtheta) d2rho_dsalin_dpress = denominator_r(i,j,k) & *( d2den_dsalin_dpress & -density_press(i,j,k)*dden_dsalinity & -density_salinity(i,j,k)*dden_dpress & -rho(i,j,k)*d2den_dsalin_dpress) d2rho_dtheta_dpress = denominator_r(i,j,k) & *( d2den_dtheta_dpress & -density_press(i,j,k)*dden_dtheta & -density_theta(i,j,k)*dden_dpress & -rho(i,j,k)*d2den_dtheta_dpress) rho_inv = Grd%tmask(i,j,k)/(rho(i,j,k) + epsln) cabbeling_param(i,j,k) = -rho_inv* & ( d2rho_dtheta2 & -2.0*d2rho_dsalin_dtheta*rhotheta_rhosalinity & +d2rho_dsalin2*rhotheta_rhosalinity**2 ) thermobaric_param(i,j,k) = -rho_inv*(d2rho_dtheta_dpress-d2rho_dsalin_dpress*rhotheta_rhosalinity) enddo enddo enddo elseif(eos_teos10) then do k=1,nk do j=jsd,jed do i=isd,ied t1 = theta(i,j,k) t2 = t1*t1 s1 = salinity(i,j,k) sp5 = sqrt(s1) p1 = press(i,j,k) - press_standard p2 = p1*p1 p1t1 = p1*t1 if(s1 > 0.0) then spm5 = 1.0/sp5 else spm5 = 0.0 endif rhotheta_rhosalinity = Grd%tmask(i,j,k)*density_theta(i,j,k)/(epsln+density_salinity(i,j,k)) dden_dtheta = v22 + t1*(two_v23 +t1*(three_v24+four_v25*t1)) & + s1*(v27 + t1*(two_v28 + t1*(three_v29 + four_v30*t1)) & + sp5*(v32+t1*(two_v33 + t1*(three_v34 + four_v35*t1)))) & + p1*((v38 + t1*(two_v39 + three_v40*t1)) & + s1*v42 & + p1*((v44 + two_v45*t1 + v46*s1) & + p1*(v48))) dden_dsalinity = & (v26 + t1*(v27 + t1*(v28 + t1*(v29 + v30*t1))) + two_v36*s1 & + 1.5*sp5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1))))) & + p1*( v41 + v42*t1 + p1*t1*v46 ) dden_dpress = (v37 + t1*(v38 + t1*(v39 + v40*t1)) & + s1*(v41 + v42*t1) & + 2.*p1*(v43 + t1*(v44 + v45*t1 + v46*s1) & + 1.5*p1*(v47 + v48*t1))) d2num_dtheta2 = (two_v03 + six_v04*t1) & + s1*(two_v07 & + sp5*( (two_v10 + six_v11*t1))) & + p1*(two_v14 & + p1*two_v19) d2num_dsalin2 = 0.75*(v08 + t1*(v09 + t1*(v10 + v11*t1)))*spm5 d2num_dtheta_dsalin = v06 + two_v07*t1 & + 1.5*sp5*(v09 + t1*(two_v10 + three_v11*t1)) & + p1*v16 d2num_dsalin_dpress = (v15 + v16*t1) & + 2.0*p1*v20 d2num_dtheta_dpress = (v13 + two_v14*t1) + s1*v16 & + 2.0*p1*(v18 + two_v19*t1 ) d2den_dtheta2 = (two_v23 +t1*(six_v24 + twelve_v25*t1)) & + s1*(two_v28 + t1*(six_v29 + twelve_v30*t1) & + sp5*((two_v33 + t1*(six_v34 + twelve_v35*t1)))) & + p1*(((two_v39 + six_v40*t1)) & + p1*two_v45 ) d2den_dsalin2 = two_v36 & + 0.75*spm5*(v31 + t1*(v32 + t1*(v33 + t1*(v34 + v35*t1)))) d2den_dsalin_dtheta = ((v27 + t1*(two_v28 + t1*(three_v29 + four_v30*t1))) & + 1.5*sp5*((v32 + t1*(two_v33 + t1*(three_v34 + four_v35*t1))))) & + p1*( v42 + p1*v46 ) d2den_dsalin_dpress = v41 + v42*t1 + p1*t1*two_v46 d2den_dtheta_dpress = (v38 + t1*(two_v39 + three_v40*t1)) & + s1*v42 & + 2.*p1*((v44 + two_v45*t1 + v46*s1) & + 1.5*p1*(v48)) d2rho_dtheta2 = denominator_r(i,j,k) & *(d2num_dtheta2-2.0*density_theta(i,j,k)*dden_dtheta-rho(i,j,k)*d2den_dtheta2) d2rho_dsalin2 = denominator_r(i,j,k) & *(d2num_dsalin2-2.0*density_salinity(i,j,k)*dden_dsalinity-rho(i,j,k)*d2den_dsalin2) d2rho_dsalin_dtheta = denominator_r(i,j,k) & *( d2num_dtheta_dsalin & -density_salinity(i,j,k)*dden_dtheta & -density_theta(i,j,k)*dden_dsalinity & -rho(i,j,k)*d2den_dsalin_dtheta) d2rho_dsalin_dpress = denominator_r(i,j,k) & *( d2den_dsalin_dpress & -density_press(i,j,k)*dden_dsalinity & -density_salinity(i,j,k)*dden_dpress & -rho(i,j,k)*d2den_dsalin_dpress) d2rho_dtheta_dpress = denominator_r(i,j,k) & *( d2den_dtheta_dpress & -density_press(i,j,k)*dden_dtheta & -density_theta(i,j,k)*dden_dpress & -rho(i,j,k)*d2den_dtheta_dpress) rho_inv = Grd%tmask(i,j,k)/(rho(i,j,k) + epsln) cabbeling_param(i,j,k) = -rho_inv* & ( d2rho_dtheta2 & -2.0*d2rho_dsalin_dtheta*rhotheta_rhosalinity & +d2rho_dsalin2*rhotheta_rhosalinity**2 ) thermobaric_param(i,j,k) = -rho_inv*(d2rho_dtheta_dpress-d2rho_dsalin_dpress*rhotheta_rhosalinity) enddo enddo enddo endif ! endif for eos choice end subroutine calc_cabbeling_thermobaricity ! NAME="calc_cabbeling_thermobaricity" !####################################################################### ! ! ! ! rho(k)-rho(k+1) for all i,j with both temperatures referenced to the ! deeper pressure depth. ! ! Of use for KPP scheme. ! ! function density_delta_z (rho_initial, salinity, theta, press) real, intent(in), dimension(isd:,jsd:,:) :: rho_initial real, intent(in), dimension(isd:,jsd:,:) :: salinity real, intent(in), dimension(isd:,jsd:,:) :: theta real, intent(in), dimension(isd:,jsd:,:) :: press real, dimension(isd:ied,jsd:jed,nk) :: density_delta_z real, dimension(isd:ied,jsd:jed,nk) :: pressure integer :: k if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_delta_z): module must be initialized') endif if(density_equal_potrho) then pressure(:,:,:) = Grd%tmask(:,:,:)*potrho_press else pressure(:,:,:) = press(:,:,:) endif do k=1,nk-1 wrk1(:,:,k) = pressure(:,:,k+1) enddo wrk1(:,:,nk) = pressure(:,:,nk) wrk2(:,:,:) = density(salinity(:,:,:), theta(:,:,:), wrk1(:,:,:)) do k=1,nk-1 density_delta_z(:,:,k) = (wrk2(:,:,k) - rho_initial(:,:,k+1))*Grd%tmask(:,:,k+1) enddo density_delta_z(:,:,nk) = (wrk2(:,:,nk) - rho_initial(:,:,nk))*Grd%tmask(:,:,nk) end function density_delta_z ! NAME="density_delta_z" !####################################################################### ! ! ! ! rho(1)-rho(k+1) for all i,j. ! ! Of use for KPP scheme. ! ! function density_delta_sfc (rho_initial, salinity, theta, press) real, intent(in), dimension(isd:,jsd:,:) :: rho_initial real, intent(in), dimension(isd:,jsd:,:) :: salinity real, intent(in), dimension(isd:,jsd:,:) :: theta real, intent(in), dimension(isd:,jsd:,:) :: press real, dimension(isd:ied,jsd:jed,nk) :: density_delta_sfc real, dimension(isd:ied,jsd:jed,nk) :: pressure integer :: k if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (density_delta_sfc): module must be initialized') endif if(density_equal_potrho) then pressure(:,:,:) = Grd%tmask(:,:,:)*potrho_press else pressure(:,:,:) = press(:,:,:) endif do k=1,nk-1 wrk1(:,:,k) = pressure(:,:,k+1) enddo wrk1(:,:,nk) = pressure(:,:,nk) wrk2(:,:,:) = density_sfc(salinity, theta, wrk1) do k=1,nk-1 density_delta_sfc(:,:,k) = (wrk2(:,:,k) - rho_initial(:,:,k+1))*Grd%tmask(:,:,k+1) enddo density_delta_sfc(:,:,nk) = (wrk2(:,:,nk) - rho_initial(:,:,nk))*Grd%tmask(:,:,nk) end function density_delta_sfc ! NAME="density_delta_sfc" !####################################################################### ! ! ! ! ! Write density and pressure_at_depth to a restart. ! ! subroutine ocean_density_end(Time, Dens, use_blobs) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens logical, intent(in) :: use_blobs integer :: stdoutunit stdoutunit=stdout() if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (ocean_density_end): module must be initialized') endif if(.not. write_a_restart) then write(stdoutunit,'(/a)') '==>Warning from ocean_density_mod (ocean_density_end): NO restart written.' call mpp_error(WARNING,'==>Warning from ocean_density_mod (ocean_density_end): NO restart written.') return endif call ocean_density_restart(Time, Dens) write(stdoutunit,'(/a)') ' From ocean_density_mod: ending density chksums' call write_timestamp(Time%model_time) call ocean_density_chksum(Time, Dens, use_blobs) module_is_initialized = .FALSE. end subroutine ocean_density_end ! NAME="ocean_density_end" !####################################################################### ! ! ! Write out restart files registered through register_restart_file ! subroutine ocean_density_restart(Time, Dens, time_stamp) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens character(len=*), intent(in), optional :: time_stamp integer :: taup1 if ( .not. module_is_initialized ) then call mpp_error(FATAL, & '==>Error in ocean_density_mod (ocean_density_end): module must be initialized') endif taup1 = Time%taup1 call reset_field_pointer(Den_restart, id_restart_rho, Dens%rho(:,:,:,taup1) ) call reset_field_pointer(Den_restart, id_restart_rho_s, Dens%rho_salinity(:,:,:,taup1) ) call save_restart(Den_restart, time_stamp) end subroutine ocean_density_restart ! NAME="ocean_density_restart" !####################################################################### ! ! ! ! Compute checksums for density. ! subroutine ocean_density_chksum(Time, Dens, use_blobs) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens logical, intent(in) :: use_blobs integer :: taup1 if (.not. module_is_initialized) then call mpp_error(FATAL,& '==>Error from ocean_density_mod (ocean_density_chksum): module not yet initialized ') endif taup1 = Time%taup1 call write_chksum_3d('rho(taup1)', Dens%rho(COMP,:,taup1)*Grd%tmask(COMP,:)) if (use_blobs) then call write_chksum_3d('rhoT(taup1)', Dens%rhoT(COMP,:)*Grd%tmask(COMP,:)) endif call write_chksum_3d('pressure_at_depth', Dens%pressure_at_depth(COMP,:)*Grd%tmask(COMP,:)) call write_chksum_3d('denominator_r', denominator_r(COMP,:)*Grd%tmask(COMP,:)) call write_chksum_3d('drhodT', Dens%drhodT(COMP,:)*Grd%tmask(COMP,:)) call write_chksum_3d('drhodS', Dens%drhodS(COMP,:)*Grd%tmask(COMP,:)) call write_chksum_3d('drhodz_zt', Dens%drhodz_zt(COMP,:)*Grd%tmask(COMP,:)) end subroutine ocean_density_chksum ! NAME="ocean_density_chksum" !####################################################################### ! ! ! ! Diagnose the buoyancy frequency, both at T-points and at ! vertical interfaces of T-cells. ! ! Author: Stephen.Griffies ! ! ! subroutine compute_buoyfreq(Time, Thickness, salinity, theta, Dens, use_blobs) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta type(ocean_density_type), intent(inout) :: Dens logical, intent(in) :: use_blobs integer :: i,j,k,kp1,kbot,m integer :: tau, taup1 real :: drhodz_prev, tmpdrhodz if (.not. module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_density (compute_buoyfreq): module needs initialization') endif tau = Time%tau taup1 = Time%taup1 wrk1(:,:,:) = 0.0 wrk2(:,:,:) = 0.0 wrk3(:,:,:) = 0.0 wrk4(:,:,:) = 0.0 if (use_blobs) then ! vertical derivative of temp and salt do k=1,nk-1 kp1 = k+1 do j=jsd,jed do i=isd,ied wrk1(i,j,k) = (theta(i,j,k)-theta(i,j,kp1))/Thickness%dzwtT(i,j,k) wrk2(i,j,k) = (salinity(i,j,k)-salinity(i,j,kp1))/Thickness%dzwtT(i,j,k) enddo enddo enddo else ! vertical derivative of temp and salt do k=1,nk-1 kp1 = k+1 do j=jsd,jed do i=isd,ied wrk1(i,j,k) = (theta(i,j,k)-theta(i,j,kp1))/Thickness%dzwt(i,j,k) wrk2(i,j,k) = (salinity(i,j,k)-salinity(i,j,kp1))/Thickness%dzwt(i,j,k) enddo enddo enddo endif ! vertical derivative of locally referenced potential density. ! vanishes at the bottom of a column. do k=1,nk-1 kp1=k+1 do j=jsd,jed do i=isd,ied wrk3(i,j,k) = Dens%drhodT(i,j,k)*wrk1(i,j,k) + Dens%drhodS(i,j,k)*wrk2(i,j,k) if(abs(wrk3(i,j,k)) < epsln_drhodz) then wrk3(i,j,k) = sign(1.0,wrk3(i,j,k))*epsln_drhodz endif wrk3(i,j,k) = Grd%tmask(i,j,kp1)*wrk3(i,j,k) enddo enddo enddo ! vertically smooth drhodz; otherwise can get noisy results if(buoyfreq_smooth_vert) then do m=1,num_121_passes do j=jsd,jed do i=isd,ied drhodz_prev = onefourth*wrk3(i,j,1) kbot=Grd%kmt(i,j) if(kbot>3) then do k=2,kbot-2 tmpdrhodz = wrk3(i,j,k) wrk3(i,j,k) = drhodz_prev + onehalf*wrk3(i,j,k) + onefourth*wrk3(i,j,k+1) drhodz_prev = onefourth*tmpdrhodz enddo endif enddo enddo enddo endif ! drhodz at bottom of T-cell do k=1,nk do j=jsd,jed do i=isd,ied Dens%drhodz_wt(i,j,k) = wrk3(i,j,k) enddo enddo enddo if(id_buoyfreq2_wt > 0) then call diagnose_3d(Time, Grd, id_buoyfreq2_wt, -grav*rho0r*wrk3(:,:,:)) endif call diagnose_3d(Time, Grd, id_drhodz_wt, wrk3(:,:,:)) ! drhodz and squared buoyancy frequency at T-cell point wrk4(:,:,:) = 0.0 k=1 do j=jsd,jed do i=isd,ied Dens%drhodz_zt(i,j,k) = Dens%drhodz_wt(i,j,k) wrk4(i,j,k) = Dens%drhodz_wt(i,j,k) enddo enddo do k=2,nk do j=jsd,jed do i=isd,ied wrk4(i,j,k) = wrk3(i,j,k-1)*(Thickness%depth_zwt(i,j,k)-Thickness%depth_zt(i,j,k)) & +wrk3(i,j,k) *(Thickness%depth_zt(i,j,k) -Thickness%depth_zwt(i,j,k-1)) wrk4(i,j,k) = wrk4(i,j,k)/(epsln+Thickness%dzt(i,j,k)) Dens%drhodz_zt(i,j,k) = wrk4(i,j,k) enddo enddo enddo if(drhodz_diag_stable) then ! enforce stable stratification do k=1,nk do j=jsd,jed do i=isd,ied if(abs(Dens%drhodz_zt(i,j,k)) < epsln_drhodz_diag .or. Dens%drhodz_zt(i,j,k) > 0.0) then Dens%drhodz_diag(i,j,k) = -epsln_drhodz_diag else Dens%drhodz_diag(i,j,k) = Dens%drhodz_zt(i,j,k) endif enddo enddo enddo else ! allow for unstable stratification do k=1,nk do j=jsd,jed do i=isd,ied if(abs(Dens%drhodz_zt(i,j,k)) < epsln_drhodz_diag) then Dens%drhodz_diag(i,j,k) = sign(1.0,Dens%drhodz_zt(i,j,k))*epsln_drhodz_diag else Dens%drhodz_diag(i,j,k) = Dens%drhodz_zt(i,j,k) endif enddo enddo enddo endif if(id_buoyfreq2_zt > 0) then call diagnose_3d(Time, Grd, id_buoyfreq2_zt, -grav*rho0r*wrk4(:,:,:)) endif call diagnose_3d(Time, Grd, id_drhodz_zt, wrk4(:,:,:)) call diagnose_3d(Time, Grd, id_drhodz_diag, Dens%drhodz_diag(:,:,:)) ! dTdz and dSdz at T-cell point wrk3(:,:,:) = 0.0 wrk4(:,:,:) = 0.0 k=1 do j=jsd,jed do i=isd,ied Dens%dTdz_zt(i,j,k) = wrk1(i,j,k) Dens%dSdz_zt(i,j,k) = wrk2(i,j,k) enddo enddo do k=2,nk do j=jsd,jed do i=isd,ied wrk3(i,j,k) = wrk1(i,j,k-1)*(Thickness%depth_zwt(i,j,k)-Thickness%depth_zt(i,j,k)) & +wrk1(i,j,k) *(Thickness%depth_zt(i,j,k) -Thickness%depth_zwt(i,j,k-1)) wrk4(i,j,k) = wrk2(i,j,k-1)*(Thickness%depth_zwt(i,j,k)-Thickness%depth_zt(i,j,k)) & +wrk2(i,j,k) *(Thickness%depth_zt(i,j,k) -Thickness%depth_zwt(i,j,k-1)) wrk3(i,j,k) = wrk3(i,j,k)/(epsln+Thickness%dzt(i,j,k)) wrk4(i,j,k) = wrk4(i,j,k)/(epsln+Thickness%dzt(i,j,k)) Dens%dTdz_zt(i,j,k) = wrk3(i,j,k) Dens%dSdz_zt(i,j,k) = wrk4(i,j,k) enddo enddo enddo end subroutine compute_buoyfreq ! NAME="compute_buoyfreq" !####################################################################### ! ! ! ! Diagnose the square of the buoyancy frequency at the bottom of ! T-cells, NOT at T-points. ! The algorithm follows that used by the private function ! compute_buoyfreq in the density module. ! ! We take the square of the buoyancy frequency as is, we do not ! smooth or force it to be positive. This allows us to search for ! instabilities. ! ! Authors: m.bates ! ! function buoyfreq2(Time, Thickness, Dens, salinity, theta, rho, use_blobs) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(in) :: rho logical, intent(in) :: use_blobs real, dimension(isd:ied,jsd:jed,nk) :: buoyfreq2 real, dimension(isd:ied,jsd:jed,nk) :: drhodT real, dimension(isd:ied,jsd:jed,nk) :: drhodS integer :: i,j,k,kp1 integer :: taup1 real :: tmp, grav_rho0r if (.not. module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_density (buoyfreq2): module needs initialization') endif taup1 = Time%taup1 grav_rho0r = grav*rho0r wrk4(:,:,:) = 0.0 buoyfreq2(:,:,:) = 0.0 !get the partial derivatives drho/dT and drho/dS call density_derivs( rho(:,:,:), salinity(:,:,:), theta(:,:,:), & Dens%pressure_at_depth(:,:,:), & Time, drhodT(:,:,:), drhodS(:,:,:), & wrk4(:,:,:)) wrk3(:,:,:) = 0.0 ! We don't need drho/dp and so can overwrite it. wrk4(:,:,:) = 0.0 ! vertical derivative of temp and salt if (use_blobs) then do k=1,nk-1 kp1 = k+1 do j=jsd,jed do i=isd,ied tmp = Grd%tmask(i,j,kp1)/Thickness%dzwtT(i,j,k) wrk3(i,j,k) = tmp*(theta(i,j,k) - theta(i,j,kp1) ) wrk4(i,j,k) = tmp*(salinity(i,j,k) - salinity(i,j,kp1)) enddo enddo enddo else do k=1,nk-1 kp1 = k+1 do j=jsd,jed do i=isd,ied tmp = Grd%tmask(i,j,kp1)/Thickness%dzwt(i,j,k) wrk3(i,j,k) = tmp*(theta(i,j,k) - theta(i,j,kp1) ) wrk4(i,j,k) = tmp*(salinity(i,j,k) - salinity(i,j,kp1)) enddo enddo enddo endif ! vertical derivative of locally referenced potential density ! at the bottom of a cell. ! N**2==0 at the bottom of a column. do k=1,nk-1 do j=jsd,jed do i=isd,ied buoyfreq2(i,j,k) = -grav_rho0r*( drhodT(i,j,k)*wrk3(i,j,k) + drhodS(i,j,k)*wrk4(i,j,k) ) enddo enddo enddo end function buoyfreq2 ! NAME="buoyfreq2" !####################################################################### ! ! ! ! Diagnose horiz derivatives of locally referenced potential density. ! For use in diagnosing balanced Ertel PV. ! ! For simplicity, do no spatial averaging. Hence, the placement ! of the derivative is at the cell face, not cell center. But to ! simplify matters, we consider the derivs to be at the cell center, ! which simplifies how we take |grad b| for the PV calculation, meaning ! there will be no averaging operations. This assumption may need to ! be revisited. ! ! Author: Stephen.Griffies ! ! ! subroutine compute_drhodxy(Time, salinity, theta, Dens) type(ocean_time_type), intent(in) :: Time real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta type(ocean_density_type), intent(inout) :: Dens integer :: i,j,k real :: tmp1,tmp2 if (.not. module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_density (compute_drhodxy): module needs initialization') endif wrk1(:,:,:) = 0.0 wrk2(:,:,:) = 0.0 ! x-derivative do k=1,nk do j=jsd,jed do i=isd,iec tmp1 = (theta(i+1,j,k) -theta(i,j,k) )*Dens%drhodT(i,j,k) tmp2 = (salinity(i+1,j,k)-salinity(i,j,k))*Dens%drhodS(i,j,k) wrk1(i,j,k) = Grd%tmask(i+1,j,k)*(tmp1+tmp2)*Grd%dxter(i,j) Dens%drhodx_zt(i,j,k) = wrk1(i,j,k) enddo enddo enddo ! y-derivative do k=1,nk do j=jsd,jec do i=isd,ied tmp1 = (theta(i,j+1,k) -theta(i,j,k) )*Dens%drhodT(i,j,k) tmp2 = (salinity(i,j+1,k)-salinity(i,j,k))*Dens%drhodS(i,j,k) wrk2(i,j,k) = Grd%tmask(i,j+1,k)*(tmp1+tmp2)*Grd%dytnr(i,j) Dens%drhody_zt(i,j,k) = wrk2(i,j,k) enddo enddo enddo if(id_drhodx_zt > 0) then call diagnose_3d(Time, Grd, id_drhodx_zt, wrk1(:,:,:)) endif if(id_drhody_zt > 0) then call diagnose_3d(Time, Grd, id_drhody_zt, wrk2(:,:,:)) endif end subroutine compute_drhodxy ! NAME="compute_drhodxy" end module ocean_density_mod