diff --git a/src/mom5/ocean_core/ocean_barotropic.F90 b/src/mom5/ocean_core/ocean_barotropic.F90 index bde1f2b604..7389db3b56 100644 --- a/src/mom5/ocean_core/ocean_barotropic.F90 +++ b/src/mom5/ocean_core/ocean_barotropic.F90 @@ -480,7 +480,8 @@ module ocean_barotropic_mod use ocean_types_mod, only: ocean_adv_vel_type, ocean_prog_tracer_type use ocean_types_mod, only: ocean_lagrangian_type use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_2d_u, diagnose_2d_en, write_chksum_2d -use ocean_workspace_mod, only: wrk1_2d, wrk2_2d, wrk3_2d, wrk4_2d, wrk1_v2d, wrk2_v2d +use ocean_workspace_mod, only: wrk1_2d, wrk2_2d, wrk3_2d, wrk4_2d, wrk1_v2d, wrk2_v2d, wrk1 +use ocean_tracer_diag_mod, only: compute_budget_mld, compute_tracer_at_mlb implicit none @@ -543,7 +544,6 @@ module ocean_barotropic_mod logical :: verbose_init=.true. ! for printouts during initialization integer :: diag_step=-1 ! for printing ascii diagnostics - ! for vertical coordinate integer :: vert_coordinate integer :: vert_coordinate_class @@ -771,6 +771,12 @@ module ocean_barotropic_mod integer :: id_eta_u =-1 integer :: id_eta_t_bar =-1 integer :: id_eta_t_tendency =-1 +integer :: id_eta_t_tendency_times_temp_in_mld =-1 +integer :: id_eta_t_tendency_times_salt_in_mld =-1 +integer :: id_eta_t_tendency_times_temp_at_mlb =-1 +integer :: id_eta_t_tendency_times_salt_at_mlb =-1 +integer :: id_s_surf_ent_temp =-1 +integer :: id_s_surf_ent_salt =-1 integer :: id_smooth_lap =-1 integer :: id_smooth_lap_diag =-1 integer :: id_smooth_bih =-1 @@ -848,6 +854,8 @@ module ocean_barotropic_mod integer :: id_grad_anompby =-1 integer :: id_eta_smoother =-1 +integer :: id_eta_smoother_times_temp_in_mld =-1 +integer :: id_eta_smoother_times_salt_in_mld =-1 integer :: id_pbot_smoother =-1 integer :: id_smooth_mask =-1 @@ -2015,6 +2023,16 @@ subroutine barotropic_diag_init(Time) Time%model_time,'surface smoother applied to eta', 'm/s', & missing_value=missing_value, range=(/-1e6,1e6/)) + id_eta_smoother_times_temp_in_mld = register_diag_field ('ocean_model', 'eta_smoother_times_temp_in_mld', & + Grd%tracer_axes(1:2), & + Time%model_time, 'tendency for eta_t over a time step times temp in mld', 'deg_C/s', & + missing_value=missing_value, range=(/-1e6,1e6/)) + + id_eta_smoother_times_salt_in_mld = register_diag_field ('ocean_model', 'eta_smoother_times_salt_in_mld', & + Grd%tracer_axes(1:2), & + Time%model_time, 'tendency for eta_t over a time step times salt in mld', 'psu kg m-3 s-1', & + missing_value=missing_value, range=(/-1e6,1e6/)) + id_pbot_smoother = register_diag_field ('ocean_model', 'pbot_smooth', Grd%tracer_axes(1:2), & Time%model_time,'bottom smoother applied to bottom pressure', 'dbar/s', & missing_value=missing_value, range=(/-1e6,1e6/)) @@ -2076,6 +2094,36 @@ subroutine barotropic_diag_init(Time) Time%model_time, 'tendency for eta_t over a time step', 'm/s', & missing_value=missing_value, range=(/-1e4,1e4/)) + id_eta_t_tendency_times_temp_in_mld = register_diag_field ('ocean_model', 'eta_t_tendency_times_temp_in_mld', & + Grd%tracer_axes(1:2), & + Time%model_time, 'tendency for eta_t over a time step times temp in mld', 'kg m-3 deg_C/s', & + missing_value=missing_value, range=(/-1e6,1e6/)) + + id_eta_t_tendency_times_salt_in_mld = register_diag_field ('ocean_model', 'eta_t_tendency_times_salt_in_mld', & + Grd%tracer_axes(1:2), & + Time%model_time, 'tendency for eta_t over a time step times salt in mld', 'psu kg m-3 s-1', & + missing_value=missing_value, range=(/-1e6,1e6/)) + + id_eta_t_tendency_times_temp_at_mlb = register_diag_field ('ocean_model', 'eta_t_tendency_times_temp_at_mlb', & + Grd%tracer_axes(1:2), & + Time%model_time, 'tendency for eta_t over a time step times temp at mlb', 'kg m-3 deg_C/s', & + missing_value=missing_value, range=(/-1e6,1e6/)) + + id_eta_t_tendency_times_salt_at_mlb = register_diag_field ('ocean_model', 'eta_t_tendency_times_salt_at_mlb', & + Grd%tracer_axes(1:2), & + Time%model_time, 'tendency for eta_t over a time step times salt at mlb', 'psu kg m-3 s-1', & + missing_value=missing_value, range=(/-1e6,1e6/)) + + id_s_surf_ent_temp = register_diag_field ('ocean_model', 's_surf_ent_temp', & + Grd%tracer_axes(1:2), & + Time%model_time, 'warming due to entrainment across s-surface at base of mixed layer', 'kg m-3 deg_C/s', & + missing_value=missing_value, range=(/-1e6,1e6/)) + + id_s_surf_ent_salt = register_diag_field ('ocean_model', 's_surf_ent_salt', & + Grd%tracer_axes(1:2), & + Time%model_time, 'salinification due to entrainment across s-surface at base of mixed layer', & + 'psu kg m-3 s-1', missing_value=missing_value, range=(/-1e6,1e6/)) + id_udrho_bt_lap = register_diag_field ('ocean_model', 'udrho_bt_lap', & Grd%vel_axes_u(1:2), Time%model_time, & 'laplacian friction to i-component udrho on barotropic time step', & @@ -2145,7 +2193,7 @@ subroutine barotropic_diag_init(Time) rho_water_r=rho0r id_pme_velocity = register_diag_field('ocean_model','pme_velocity', Grd%tracer_axes(1:2), & - Time%model_time, 'net (precip-evap)(kg/(m2*sec) into ocean, divided by 1035kg/m3', 'm/sec',& + Time%model_time, 'net (precip-evap)(kg/(m2*sec)) into ocean, divided by 1035kg/m3', 'm/sec',& missing_value=missing_value,range=(/-1e6,1e6/)) id_conv_ud_pred_bt1 = register_diag_field ('ocean_model', 'conv_ud_pred_bt1', Grd%tracer_axes(1:2), & @@ -2406,12 +2454,13 @@ end subroutine eta_and_pbot_update ! ! ! -subroutine eta_and_pbot_diagnose (Time, Dens, Thickness, patm, pme, river, Ext_mode, & +subroutine eta_and_pbot_diagnose (Time, Dens, Thickness, T_prog, patm, pme, river, Ext_mode, & L_system, use_blobs) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens type(ocean_thickness_type), intent(inout) :: Thickness + type(ocean_prog_tracer_type), intent(inout) :: T_prog(:) type(ocean_external_mode_type), intent(inout) :: Ext_mode type(ocean_lagrangian_type), intent(in) :: L_system real, dimension(isd:,jsd:), intent(in) :: patm @@ -2419,6 +2468,15 @@ subroutine eta_and_pbot_diagnose (Time, Dens, Thickness, patm, pme, river, Ext_m real, dimension(isd:,jsd:), intent(in) :: river logical, intent(in) :: use_blobs + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d + real, dimension(isd:ied,jsd:jed) :: tracer_in_mld + real, dimension(isd:ied,jsd:jed) :: mld + real, dimension(isd:ied,jsd:jed) :: eta_t_tendency + real, dimension(isd:ied,jsd:jed,1:nk) :: tracer_ext + + integer :: num_prog_tracers, index_temp, index_salt, n + integer :: tau, taup1 integer :: i,j,k,km1 real :: eta_global @@ -2426,8 +2484,14 @@ subroutine eta_and_pbot_diagnose (Time, Dens, Thickness, patm, pme, river, Ext_m tau = Time%tau taup1 = Time%taup1 - ! diagnose bottom pressure - if(vert_coordinate_class==DEPTH_BASED) then + 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 + enddo + + ! diagnose bottom pressure + if(vert_coordinate_class==DEPTH_BASED) then Ext_mode%pbot_t(:,:,taup1) = 0.0 if (use_blobs) then do k=1,nk @@ -2638,15 +2702,114 @@ subroutine eta_and_pbot_diagnose (Time, Dens, Thickness, patm, pme, river, Ext_m call diagnose_2d(Time, Grd, id_patm_for_sea_lev, Ext_mode%patm_for_sea_lev(:,:)) call diagnose_2d(Time, Grd, id_sea_lev_for_coupler, Thickness%sea_lev(:,:)) - if(id_eta_t_tendency > 0) then - wrk1_2d(:,:)= 0.0 + if(id_eta_t_tendency > 0 .or. id_eta_t_tendency_times_temp_in_mld > 0 .or. id_eta_t_tendency_times_salt_in_mld > 0 & + .or. id_eta_t_tendency_times_temp_at_mlb > 0 .or. id_eta_t_tendency_times_salt_at_mlb > 0 & + .or. id_s_surf_ent_temp > 0 .or. id_s_surf_ent_salt > 0 ) then + wrk1_2d(:,:)= 0.0 do j=jsc,jec do i=isc,iec wrk1_2d(i,j) = (Ext_mode%eta_t(i,j,taup1)-Ext_mode%eta_t(i,j,tau))*dtimer enddo enddo - call diagnose_2d(Time, Grd, id_eta_t_tendency, wrk1_2d(:,:)) - endif + eta_t_tendency(:,:) = wrk1_2d(:,:) + + if (id_eta_t_tendency > 0) then + call diagnose_2d(Time, Grd, id_eta_t_tendency, wrk1_2d(:,:)) + endif + + if ( id_eta_t_tendency_times_temp_in_mld > 0) then + wrk1(:,:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1(i,j,k) = T_prog(index_temp)%field(i,j,k,tau)*Thickness%rho_dzt(i,j,k,tau) + enddo + enddo + enddo + tracer_ext(:,:,:) = wrk1(:,:,:) + tracer_in_mld(:,:) = 0.0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, tracer_ext(:,:,:), tracer_in_mld(:,:)) + + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = eta_t_tendency(:,:)*tracer_in_mld(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_eta_t_tendency_times_temp_in_mld, tendency_in_mld(:,:)) + endif + + if ( id_eta_t_tendency_times_temp_at_mlb > 0 .or. id_s_surf_ent_temp > 0 ) then + tracer_ext(:,:,:) = T_prog(index_temp)%field(:,:,:,tau) + tracer_in_mld(:,:) = 0.0 + mld(:,:) = 0.0 + call compute_tracer_at_mlb(Time, Thickness, Dens, T_prog, tracer_ext(:,:,:), tracer_in_mld(:,:), mld(:,:)) + + tendency_3d(:,:,:) = 0.0 + + if ( id_eta_t_tendency_times_temp_at_mlb > 0 ) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,1) = eta_t_tendency(:,:)*tracer_in_mld(:,:)*rho0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_eta_t_tendency_times_temp_at_mlb, tendency_in_mld(:,:)) + endif + + if ( id_s_surf_ent_temp > 0 ) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,1) = Grd%tmask(:,:,1)*rho0*eta_t_tendency(:,:)*tracer_in_mld(:,:)*(1.0 - & + (mld(:,:)/(Grd%ht(:,:) + Ext_mode%eta_t(:,:,tau) + epsln)) & + ) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_s_surf_ent_temp, tendency_in_mld(:,:)) + endif + + endif + + if ( id_eta_t_tendency_times_salt_in_mld > 0) then + wrk1(:,:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1(i,j,k) = T_prog(index_salt)%field(i,j,k,tau)*Thickness%rho_dzt(i,j,k,tau) + enddo + enddo + enddo + tracer_ext(:,:,:) = wrk1(:,:,:) + tracer_in_mld(:,:) = 0.0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, tracer_ext(:,:,:), tracer_in_mld(:,:)) + + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = eta_t_tendency(:,:)*tracer_in_mld(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_eta_t_tendency_times_salt_in_mld, tendency_in_mld(:,:)) + endif + + if ( id_eta_t_tendency_times_salt_at_mlb > 0 .or. id_s_surf_ent_salt > 0 ) then + tracer_ext(:,:,:) = T_prog(index_salt)%field(:,:,:,tau) + tracer_in_mld(:,:) = 0.0 + mld(:,:) = 0.0 + call compute_tracer_at_mlb(Time, Thickness, Dens, T_prog, tracer_ext(:,:,:), tracer_in_mld(:,:), mld(:,:)) + + tendency_3d(:,:,:) = 0.0 + + if ( id_eta_t_tendency_times_salt_at_mlb > 0 ) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,1) = eta_t_tendency(:,:)*tracer_in_mld(:,:)*rho0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_eta_t_tendency_times_salt_at_mlb, tendency_in_mld(:,:)) + endif + + if ( id_s_surf_ent_salt > 0 ) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,1) = Grd%tmask(:,:,1)*rho0*eta_t_tendency(:,:)*tracer_in_mld(:,:)*(1.0 - & + (mld(:,:)/(Grd%ht(:,:) + Ext_mode%eta_t(:,:,tau) + epsln)) & + ) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_s_surf_ent_salt, tendency_in_mld(:,:)) + endif + + endif + + endif end subroutine eta_and_pbot_diagnose @@ -4748,32 +4911,43 @@ end subroutine eta_smooth_diagnosed ! ! ! -! Compute tendency for smoothing eta and tracer. +! Compute tendency for smoothing eta and tracer. ! -! Use either a laplacian or a biharmonic smoothing operator. +! Use either a laplacian or a biharmonic smoothing operator. ! ! This smoothing option is most useful for the Bgrid, which has -! a checkerboard null mode when discretizing gravity waves. The +! a checkerboard null mode when discretizing gravity waves. The ! Cgrid has no gravity wave noise so may not need this smoother. ! -! Recommend against using the biharmonic, since it is not a -! positive definite operator and so can lead to extrema. -! Biharmonic is retained for legacy purposes. +! Recommend against using the biharmonic, since it is not a +! positive definite operator and so can lead to extrema. +! Biharmonic is retained for legacy purposes. ! ! ! -subroutine ocean_eta_smooth(Time, Thickness, Ext_mode, T_prog) +subroutine ocean_eta_smooth(Time, Thickness, Dens, Ext_mode, T_prog) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(inout) :: Thickness + type(ocean_density_type), intent(in) :: Dens type(ocean_external_mode_type), intent(inout) :: Ext_mode type(ocean_prog_tracer_type), intent(inout) :: T_prog(:) real, dimension(isd:ied,jsd:jed) :: tmp - integer :: i, j, n, nprog, taum1 + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d + real, dimension(isd:ied,jsd:jed) :: tracer_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tracer_ext + + integer :: i, j, n, k, nprog, taum1, tau + integer :: index_temp, index_salt real :: eta_min - nprog = size(T_prog(:)) + nprog = size(T_prog(:)) + do n=1,nprog + if (T_prog(n)%name == 'temp') index_temp = n + if (T_prog(n)%name == 'salt') index_salt = n + enddo ! initialise some fields for later use with OBC and return if(.not. smooth_eta_t_laplacian .and. .not. smooth_eta_t_biharmonic) then @@ -4787,6 +4961,7 @@ subroutine ocean_eta_smooth(Time, Thickness, Ext_mode, T_prog) endif taum1 = Time%taum1 + tau = Time%tau smooth_mask(:,:) = 0.0 etastar(:,:) = 0.0 tmp(:,:) = 0.0 @@ -4848,10 +5023,50 @@ subroutine ocean_eta_smooth(Time, Thickness, Ext_mode, T_prog) if (id_eta_smoother > 0) call diagnose_2d(Time, Grd, id_eta_smoother, tmp(:,:)*rho0r) + if (id_eta_smoother_times_temp_in_mld > 0) then + wrk1(:,:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1(i,j,k) = T_prog(index_temp)%field(i,j,k,tau)*Thickness%rho_dzt(i,j,k,tau) + enddo + enddo + enddo + tracer_ext(:,:,:) = wrk1(:,:,:) + tracer_in_mld(:,:) = 0.0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, tracer_ext(:,:,:), tracer_in_mld(:,:)) + + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = Ext_mode%eta_smooth(:,:)*rho0r*tracer_in_mld(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_eta_smoother_times_temp_in_mld, tendency_in_mld(:,:)) + endif + + if (id_eta_smoother_times_salt_in_mld > 0) then + wrk1(:,:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1(i,j,k) = T_prog(index_salt)%field(i,j,k,tau)*Thickness%rho_dzt(i,j,k,tau) + enddo + enddo + enddo + tracer_ext(:,:,:) = wrk1(:,:,:) + tracer_in_mld(:,:) = 0.0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, tracer_ext(:,:,:), tracer_in_mld(:,:)) + + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = Ext_mode%eta_smooth(:,:)*rho0r*tracer_in_mld(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_eta_smoother_times_salt_in_mld, tendency_in_mld(:,:)) + endif + ! T_prog%eta_smooth has dimensions tracer concentration * (kg/m^3)*(m/s). - ! note that tracer filter is zero when eta_t is zero, as we wish since in + ! note that tracer filter is zero when eta_t is zero, as we wish since in ! this case there is no need to apply a filter. - if(smooth_eta_t_laplacian) then + if(smooth_eta_t_laplacian) then do n=1,nprog do j=jsd,jed do i=isd,ied diff --git a/src/mom5/ocean_core/ocean_model.F90 b/src/mom5/ocean_core/ocean_model.F90 index 176773177e..29308dbd54 100644 --- a/src/mom5/ocean_core/ocean_model.F90 +++ b/src/mom5/ocean_core/ocean_model.F90 @@ -1604,7 +1604,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, & ! compute "flux adjustments" (e.g., surface tracer restoring, flux correction) call mpp_clock_begin(id_flux_adjust) - call flux_adjust(Time, T_diag(1:num_diag_tracers), Dens, Ext_mode, & + call flux_adjust(Time, T_diag(1:num_diag_tracers), Dens, Thickness, Ext_mode, & T_prog(1:num_prog_tracers), Velocity, river, melt, pme, runoff) call mpp_clock_end(id_flux_adjust) @@ -1718,7 +1718,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, & ! smoother is not needed for C-grid. if(vert_coordinate_class==DEPTH_BASED) then call mpp_clock_begin(id_surface_smooth) - call ocean_eta_smooth(Time, Thickness, Ext_mode, T_prog(1:num_prog_tracers)) + call ocean_eta_smooth(Time, Thickness, Dens, Ext_mode, T_prog(1:num_prog_tracers)) call mpp_clock_end(id_surface_smooth) endif @@ -1879,7 +1879,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, & ! diagnose time=taup1 ocean free surface height or bottom pressure. ! also diagnose geodepth_zt and geodepth_zwt. call mpp_clock_begin(id_eta_and_pbot_diagnose) - call eta_and_pbot_diagnose(Time, Dens, Thickness, patm, pme, river, Ext_mode, Lagrangian_system, use_blobs) + call eta_and_pbot_diagnose(Time, Dens, Thickness, T_prog, patm, pme, river, Ext_mode, Lagrangian_system, use_blobs) call mpp_clock_end(id_eta_and_pbot_diagnose) ! diagnose the geodepth of new blobs and the depth of old blobs diff --git a/src/mom5/ocean_core/ocean_sbc.F90 b/src/mom5/ocean_core/ocean_sbc.F90 index 4d40cecd3e..600a2f5e50 100644 --- a/src/mom5/ocean_core/ocean_sbc.F90 +++ b/src/mom5/ocean_core/ocean_sbc.F90 @@ -533,6 +533,7 @@ module ocean_sbc_mod use ocean_workspace_mod, only: wrk1_2d, wrk2_2d, wrk3_2d, wrk1 use ocean_util_mod, only: diagnose_2d, diagnose_2d_u, diagnose_3d_u, diagnose_sum use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld #if defined(CSIRO_BGC) use csiro_bgc_mod, only: csiro_bgc_virtual_fluxes, do_csiro_bgc,ind_no3,ind_phy @@ -592,6 +593,8 @@ module ocean_sbc_mod integer, allocatable, dimension(:) :: id_stf_calving integer, allocatable, dimension(:) :: id_stf_pme integer, allocatable, dimension(:) :: id_stf_pme_on_nrho +integer, allocatable, dimension(:) :: id_stf_pme_in_mld +integer, allocatable, dimension(:) :: id_pme_river_times_tracer_in_mld integer, allocatable, dimension(:) :: id_stf_prec integer, allocatable, dimension(:) :: id_stf_evap integer, allocatable, dimension(:) :: id_trunoff @@ -625,6 +628,7 @@ module ocean_sbc_mod integer :: id_wavlen =-1 integer :: id_net_sfc_heating =-1 +integer :: id_net_sfc_heating_in_mld =-1 integer :: id_total_net_sfc_heating =-1 integer :: id_net_sfc_workq =-1 @@ -641,13 +645,17 @@ module ocean_sbc_mod integer :: id_ekman_we =-1 integer :: id_ekman_heat =-1 integer :: id_swflx =-1 +integer :: id_swflx_in_mld =-1 integer :: id_swflx_vis =-1 integer :: id_lw_heat =-1 integer :: id_sens_heat =-1 +integer :: id_lw_heat_in_mld =-1 +integer :: id_sens_heat_in_mld =-1 integer :: id_fprec_melt_heat =-1 integer :: id_calving_melt_heat =-1 integer :: id_evap_heat =-1 +integer :: id_evap_heat_in_mld =-1 integer :: id_fprec =-1 integer :: id_lprec =-1 @@ -664,6 +672,7 @@ module ocean_sbc_mod integer :: id_evap =-1 integer :: id_pme_sbc =-1 integer :: id_pme_river =-1 +integer :: id_pme_river_in_mld =-1 integer :: id_pme_restore =-1 integer :: id_pme_eta_restore=-1 integer :: id_pme_correct =-1 @@ -1667,6 +1676,8 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) allocate( id_stf_calving (num_prog_tracers) ) allocate( id_stf_pme (num_prog_tracers) ) allocate( id_stf_pme_on_nrho (num_prog_tracers) ) + allocate( id_stf_pme_in_mld (num_prog_tracers) ) + allocate( id_pme_river_times_tracer_in_mld (num_prog_tracers) ) allocate( id_stf_prec (num_prog_tracers) ) allocate( id_stf_evap (num_prog_tracers) ) allocate( id_trunoff (num_prog_tracers) ) @@ -1691,6 +1702,8 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) id_stf_calving (:) = -1 id_stf_pme (:) = -1 id_stf_pme_on_nrho (:) = -1 + id_stf_pme_in_mld (:) = -1 + id_pme_river_times_tracer_in_mld (:) = -1 id_stf_prec (:) = -1 id_stf_evap (:) = -1 id_trunoff (:) = -1 @@ -1866,6 +1879,11 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) '(kg/m^3)*(m/sec)', missing_value=missing_value,range=(/-1e6,1e6/), & standard_name='water_flux_into_sea_water') + id_pme_river_in_mld= register_diag_field('ocean_model','pme_river_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'mass flux of precip-evap+river via sbc (liquid, frozen, evaporation) averaged in mixed layer',& + '(kg/m^3)*(m/sec)/m', missing_value=missing_value,range=(/-1e6,1e6/), & + standard_name='water_flux_into_sea_water') + id_net_sfc_workemp = register_diag_field('ocean_model','net_sfc_workEmP', & Grd%tracer_axes(1:2), & Time%model_time, 'pme_river*g*beta2*So/rho0, beta uses pot_rho rather than neut', & @@ -1952,6 +1970,22 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) missing_value=missing_value,range=(/-1.e10,1.e10/), & standard_name='surface_downward_sensible_heat_flux') + id_swflx_in_mld = register_diag_field('ocean_model','swflx_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'shortwave flux into ocean (>0 heats ocean) averaged in mixed layer', 'W/m^3', & + missing_value=missing_value,range=(/-1.e10,1.e10/)) + + id_evap_heat_in_mld = register_diag_field('ocean_model','evap_heat_in_mld', Grd%tracer_axes(1:2),& + Time%model_time, 'latent heat flux into ocean (<0 cools ocean) averaged in mixed layer', 'W/m^3', & + missing_value=missing_value,range=(/-1e10,1e10/)) + + id_lw_heat_in_mld = register_diag_field('ocean_model','lw_heat_in_mld', Grd%tracer_axes(1:2),& + Time%model_time, 'longwave flux into ocean (<0 cools ocean) averaged in mixed layer', 'W/m^3' , & + missing_value=missing_value,range=(/-1.e10,1.e10/)) + + id_sens_heat_in_mld = register_diag_field('ocean_model','sens_heat_in_mld', Grd%tracer_axes(1:2),& + Time%model_time, 'sensible heat into ocean (<0 cools ocean) averaged in mixed layer', 'W/m^3' , & + missing_value=missing_value,range=(/-1.e10,1.e10/)) + id_fprec_melt_heat = register_diag_field('ocean_model','fprec_melt_heat', Grd%tracer_axes(1:2),& Time%model_time, 'heat flux to melt frozen precip (<0 cools ocean)', 'W/m^2' , & missing_value=missing_value,range=(/-1.e10,1.e10/), & @@ -2396,6 +2430,11 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) Time%model_time, 'surface ocean heat flux coming through coupler and mass transfer',& 'Watts/m^2' , & missing_value=missing_value,range=(/-1.e4,1.e4/)) + id_net_sfc_heating_in_mld = register_diag_field('ocean_model','net_sfc_heating_in_mld', & + Grd%tracer_axes(1:2), & + Time%model_time, 'surface ocean heat flux coming through coupler and mass transfer averaged in mixed layer',& + 'Watts/m^3' , & + missing_value=missing_value,range=(/-1.e4,1.e4/)) id_net_sfc_workq = register_diag_field('ocean_model','net_sfc_workq', & Grd%tracer_axes(1:2), & Time%model_time, 'net_sfc_heat*g*alpha2/Cp/rho0, alpha uses pot_rho rather than neut', & @@ -2421,6 +2460,14 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) Dens%neutralrho_axes(1:3), & Time%model_time, 'heat flux (relative to 0C) from pme transfer of water across ocean surface binned to neutral density', & 'Watts/m^2' , missing_value=missing_value,range=(/-1.e20,1.e20/)) + id_stf_pme_in_mld(n) = register_diag_field('ocean_model','sfc_hflux_pme_in_mld', & + Grd%tracer_axes(1:2), & + Time%model_time, 'heat flux (relative to 0C) from pme transfer of water across ocean surface averaged in mixed layer', & + 'Watts/m^3' , missing_value=missing_value,range=(/-1.e4,1.e4/)) + id_pme_river_times_tracer_in_mld(n) = register_diag_field('ocean_model','pme_river_times_temp_in_mld', & + Grd%tracer_axes(1:2), & + Time%model_time, 'mass flux of precip-evap+river via sbc (liquid, frozen, evaporation) averaged in mixed layer times mixed layer temperature',& + '(kg/m^3)*(m/sec) '//trim(T_prog(n)%units)//' / m', missing_value=missing_value,range=(/-1e6,1e6/)) id_stf_prec(n) = register_diag_field('ocean_model','sfc_hflux_from_water_prec', & Grd%tracer_axes(1:2), & Time%model_time, 'heat flux from precip transfer of water across ocean surface',& @@ -2542,6 +2589,21 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) Time%model_time, trim(name), & 'kg/(m^2*sec)' , & missing_value=missing_value,range=(/-1.e20,1.e20/)) + name = 'sfc_'//trim(T_prog(n)%name)//'_flux_pme_in_mld' + id_stf_pme_in_mld(n) = register_diag_field('ocean_model',& + trim(name), & + Grd%tracer_axes(1:2), & + Time%model_time, trim(name), & + 'kg/(m^3*sec)' , & + missing_value=missing_value,range=(/-1.e4,1.e4/)) + name = 'pme_river_times_'//trim(T_prog(n)%name)//'_in_mld' + id_pme_river_times_tracer_in_mld(n) = register_diag_field('ocean_model',& + trim(name), & + Grd%tracer_axes(1:2), & + Time%model_time, & + 'mass flx of precip-evap+river via sbc (liquid, frozen, evaporation) averaged in mixed layer times mixed layer salinity',& + '(kg/m^3)*(m/sec) '//trim(T_prog(n)%units)//' / m',& + missing_value=missing_value,range=(/-1e6,1e6/)) name = 'sfc_'//trim(T_prog(n)%name)//'_flux_prec' id_stf_prec(n) = register_diag_field('ocean_model',& trim(name), & @@ -2669,6 +2731,13 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) Time%model_time, trim(name), & 'kg/(m^2*sec)' , & missing_value=missing_value,range=(/-1.e4,1.e4/)) + name = 'sfc_'//trim(T_prog(n)%name)//'_flux_pme_in_mld' + id_stf_pme_in_mld(n) = register_diag_field('ocean_model',& + trim(name), & + Grd%tracer_axes(1:2), & + Time%model_time, trim(name), & + 'kg/(m^3*sec)' , & + missing_value=missing_value,range=(/-1.e4,1.e4/)) name = 'sfc_'//trim(T_prog(n)%name)//'_flux_pme_on_nrho' id_stf_pme_on_nrho(n) = register_diag_field('ocean_model',& trim(name), & @@ -4451,22 +4520,22 @@ end subroutine get_ocean_sbc ! ! ! -! Subroutine to compute the surface fluxes derived from a -! restoring condition and/or correction from an input file. +! Subroutine to compute the surface fluxes derived from a +! restoring condition and/or correction from an input file. ! -! We use a convention whereby a positive -! flux enters the ocean: (+) down convention. +! We use a convention whereby a positive +! flux enters the ocean: (+) down convention. ! ! When restoring salinity, one may choose to convert this ! flux to an implied water flux, or keep it a salt flux. ! Converting to a water flux will alter the sea level, and -! so alter the concentration of other tracers. -! The default is to keep it as a salt flux. +! so alter the concentration of other tracers. +! The default is to keep it as a salt flux. ! ! ! -subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, melt, pme, runoff) +subroutine flux_adjust(Time, T_diag, Dens, Thickness, Ext_mode, T_prog, Velocity, river, melt, pme, runoff) #if defined(ACCESS_CM) || defined(ACCESS_OM) use auscom_ice_parameters_mod, only : use_ioaice, aice_cutoff @@ -4476,6 +4545,7 @@ subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, me type(ocean_time_type), intent(in) :: Time type(ocean_diag_tracer_type), intent(in) :: T_diag(:) type(ocean_density_type), intent(in) :: Dens + type(ocean_thickness_type), intent(in) :: Thickness type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_prog_tracer_type), intent(inout) :: T_prog(:) type(ocean_velocity_type), intent(inout) :: Velocity @@ -4502,6 +4572,11 @@ subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, me logical :: ice_present real :: active_cells, smftu, smftv + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d + real, dimension(isd:ied,jsd:jed) :: tracer_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tracer + #if defined(ACCESS_CM) ! Changed in CM2. Make parameter to isolate change real, parameter :: ice_detection_parameter = 0.054 @@ -4841,6 +4916,55 @@ subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, me call diagnose_2d(Time, Grd, id_stf_pme(index_temp), & pme(:,:)*T_prog(index_temp)%tpme(:,:)*T_prog(index_temp)%conversion) endif + ! heat input from net pme relative to 0 degrees C (W/m2) averaged in mixed layer + if (id_stf_pme_in_mld(index_temp) > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = pme(:,:)*T_prog(index_temp)%tpme(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_stf_pme_in_mld(index_temp), tendency_in_mld(:,:)*T_prog(index_temp)%conversion) + endif + ! pme_river times mixed layer temperature in mixed layer + if (id_pme_river_times_tracer_in_mld(index_temp) > 0) then + wrk1(:,:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1(i,j,k) = T_prog(index_temp)%field(i,j,k,tau)*Thickness%rho_dzt(i,j,k,tau) + enddo + enddo + enddo + tracer(:,:,:) = wrk1(:,:,:) + tracer_in_mld(:,:) = 0.0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, tracer(:,:,:), tracer_in_mld(:,:)) + + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = pme(:,:) + river(:,:) + tendency_3d(:,:,1) = tendency_3d(:,:,1)*tracer_in_mld(:,:)*rho0r + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_pme_river_times_tracer_in_mld(index_temp), tendency_in_mld(:,:)) + endif + ! pme_river times mixed layer salinity in mixed layer + if (id_pme_river_times_tracer_in_mld(index_salt) > 0) then + wrk1(:,:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1(i,j,k) = T_prog(index_salt)%field(i,j,k,tau)*Thickness%rho_dzt(i,j,k,tau) + enddo + enddo + enddo + tracer(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tracer(:,:,:), tracer_in_mld(:,:)) + + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = pme(:,:) + river(:,:) + tendency_3d(:,:,1) = tendency_3d(:,:,1)*tracer_in_mld(:,:)*rho0r + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_pme_river_times_tracer_in_mld(index_salt), tendency_in_mld(:,:)) + endif ! heat input from net pme relative to 0 degrees C (W/m2) binned to ! neutral density if (id_stf_pme_on_nrho(index_temp) > 0) then @@ -5271,8 +5395,10 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou real, dimension(isd:,jsd:), intent(in) :: swflx_vis real, dimension(isd:ied,jsd:jed) :: tmp_flux + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d - integer :: i,j,k + integer :: i,j,k integer :: ii, jj integer :: tau real :: total_stuff @@ -5616,7 +5742,7 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou ! the frozen precip is typically best approximated to be at 0C, rather than SST. But ! we diagnose the contribution in this manner in order to agree with the prognostic model ! methods. - if(id_net_sfc_heating > 0) then + if(id_net_sfc_heating > 0 .or. id_net_sfc_heating_in_mld > 0) then wrk1_2d(:,:) = 0.0 do j=jsc,jec do i=isc,iec @@ -5633,10 +5759,19 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou )*T_prog(index_temp)%tpme(i,j)) enddo enddo - call diagnose_2d(Time, Grd, id_net_sfc_heating, wrk1_2d(:,:)) + if (id_net_sfc_heating > 0) then + call diagnose_2d(Time, Grd, id_net_sfc_heating, wrk1_2d(:,:)) + endif + if (id_net_sfc_heating_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = wrk1_2d(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_net_sfc_heating_in_mld, tendency_in_mld(:,:)) + endif endif - - ! net_surface_heating*g*alphasfc2/rho0/Cp + + ! net_surface_heating*g*alphasfc2/rho0/Cp if(id_net_sfc_workq > 0) then wrk1_2d(:,:) = 0.0 do j=jsc,jec @@ -5681,9 +5816,17 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou ! shortwave flux (W/m2) call diagnose_2d(Time, Grd, id_swflx, swflx(:,:)) + if (id_swflx_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = swflx(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_swflx_in_mld, tendency_in_mld(:,:)) + endif + ! total shortwave heat transport (Watts) call diagnose_sum(Time, Grd, Dom, id_total_ocean_swflx, swflx, 1e-15) - ! swflx impacts on water mass transformation in neutral density classes + ! swflx impacts on water mass transformation in neutral density classes if(id_tform_rho_pbl_sw_on_nrho > 0) then wrk1(:,:,:) = 0.0 k=1 @@ -5703,7 +5846,7 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou ! evaporative heat flux (W/m2) (<0 cools ocean) - if (id_evap_heat > 0) then + if (id_evap_heat > 0 .or. id_evap_heat_in_mld > 0) then tmp_flux=0.0 do j=jsc_bnd,jec_bnd do i=isc_bnd,iec_bnd @@ -5712,9 +5855,18 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou tmp_flux(ii,jj) = -latent_heat_vapor(ii,jj)*Ice_ocean_boundary%q_flux(i,j) enddo enddo - call diagnose_2d(Time, Grd, id_evap_heat, tmp_flux(:,:)) + if (id_evap_heat > 0) then + call diagnose_2d(Time, Grd, id_evap_heat, tmp_flux(:,:)) + endif + if (id_evap_heat_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = tmp_flux(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_evap_heat_in_mld, tendency_in_mld(:,:)) + endif endif - ! total evaporative heating (Watts) + ! total evaporative heating (Watts) if (id_total_ocean_evap_heat > 0) then tmp_flux=0.0 do j=jsc_bnd,jec_bnd @@ -5744,9 +5896,17 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou ! longwave heat flux (W/m2) call diagnose_2d(Time, Grd, id_lw_heat, longwave(:,:)) - ! total longwave heating (Watts) + if (id_lw_heat_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = longwave(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_lw_heat_in_mld, tendency_in_mld(:,:)) + endif + + ! total longwave heating (Watts) call diagnose_sum(Time, Grd, Dom, id_total_ocean_lw_heat, longwave, 1e-15) - ! longwave impacts on water mass transformation in neutral density classes + ! longwave impacts on water mass transformation in neutral density classes if(id_tform_rho_pbl_lw_on_nrho > 0) then wrk1(:,:,:) = 0.0 k=1 @@ -5794,9 +5954,17 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou ! sensible heat flux (W/m2) call diagnose_2d(Time, Grd, id_sens_heat, sensible(:,:)) - ! total sensible heat transport (Watts) + if (id_sens_heat_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = sensible(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_sens_heat_in_mld, tendency_in_mld(:,:)) + endif + + ! total sensible heat transport (Watts) call diagnose_sum(Time, Grd, Dom, id_total_ocean_sens_heat, sensible, 1e-15) - ! sensible heat impacts on water mass transformation in neutral density classes + ! sensible heat impacts on water mass transformation in neutral density classes if(id_tform_rho_pbl_sens_on_nrho > 0) then wrk1(:,:,:) = 0.0 k=1 @@ -5941,12 +6109,21 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou !--------mass flux related diagnostics----------------------------- ! - ! total mass flux per area from pme and river (kg/(m2*sec)) + ! total mass flux per area from pme and river (kg/(m2*sec)) if (id_pme_river > 0) then call diagnose_2d(Time, Grd, id_pme_river, pme(:,:) + river(:,:)) endif - - !ape work from E-P+R: pme_river*g*betasfc2*So/rho0 + + ! total mass flux per area from pme and river (kg/(m2*sec)) + if (id_pme_river_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = pme(:,:) + river(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_pme_river_in_mld, tendency_in_mld(:,:)) + endif + + !ape work from E-P+R: pme_river*g*betasfc2*So/rho0 if(id_net_sfc_workemp > 0) then call diagnose_2d(Time, Grd, id_net_sfc_workemp, & grav*salinity_ref/rho0*betasfc2(:,:)*(pme(:,:) + river(:,:)) ) diff --git a/src/mom5/ocean_diag/ocean_tracer_diag.F90 b/src/mom5/ocean_diag/ocean_tracer_diag.F90 index a9846d9d1b..0b139ae9f1 100644 --- a/src/mom5/ocean_diag/ocean_tracer_diag.F90 +++ b/src/mom5/ocean_diag/ocean_tracer_diag.F90 @@ -351,6 +351,8 @@ module ocean_tracer_diag_mod public calc_potrho_mixed_layer public send_tracer_variance public diagnose_eta_tend_3dflux +public compute_budget_mld +public compute_tracer_at_mlb private compute_subduction private compute_tracer_mld @@ -694,8 +696,6 @@ subroutine ocean_tracer_diag_init(Grid, Domain, Time, Time_steps, Thickness, T_p 'Vertically integrated tracer [sum(rho_dzt*tracer)] in mixed layer for '//trim(T_prog(n)%name),& 'kg/m^2 * '//trim(T_prog(n)%units), missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_tracer_mld(n) > 0) compute_tracer_mld_diags = .true. - - enddo id_kappa_sort = register_diag_field ('ocean_model','kappa_sort', & @@ -1072,7 +1072,7 @@ subroutine calc_mixed_layer_depth(Thickness, salinity, theta, rho, pressure, & real, dimension(isd:,jsd:), intent(inout) :: hmxl logical, optional, intent(in) :: smooth_mld_input - real, parameter :: epsln=1.0e-20 ! for divisions + real, parameter :: epsln=1.0e-20 ! for divisions integer :: i, j, k, km1, kb logical :: smooth_mld_routine @@ -1672,7 +1672,6 @@ subroutine compute_tracer_mld(Time, Thickness, Dens, T_prog, salinity, theta) integer :: i,j,k,kp1,n integer :: taup1 real, dimension(isd:ied,jsd:jed) :: mld - if (.not.module_is_initialized) then call mpp_error(FATAL, & @@ -1755,7 +1754,214 @@ subroutine compute_tracer_mld(Time, Thickness, Dens, T_prog, salinity, theta) end subroutine compute_tracer_mld ! NAME="compute_tracer_mld" +!####################################################################### +! +! +! +! +! Compute the MLD-average of a 3D tracer tendency field. +! +! A per-layer weight wrk1(i,j,k) in [0,1] is built so that the weighted +! sum over k equals the depth-integrated tendency within the mixed layer. +! Three cases handle where the MLD falls relative to the layer interfaces +! (depth_zwt): +! 1. MLD shallower than bottom of layer 1: partial weight for k=1 only. +! 2. MLD within layer k (k>1): full weight for layers above, partial +! weight for the straddling layer (fraction of dzt spanned by MLD). +! 3. MLD deeper than all model levels: full weight everywhere. +! The depth-weighted column sum is then divided by the MLD to give a +! depth-averaged tendency. +! +! Ryan.Holmes, September 2024 +! +! +subroutine compute_budget_mld(Time, Thickness, Dens, T_prog, tendency, tendency_2d) + + type(ocean_time_type), intent(in) :: Time + type(ocean_thickness_type), intent(in) :: Thickness + type(ocean_density_type), intent(in) :: Dens + type(ocean_prog_tracer_type), intent(in) :: T_prog(:) + real, dimension(isd:,jsd:,:), intent(in) :: tendency ! 3D tendency field (any units/m) + real, dimension(isd:,jsd:), intent(inout) :: tendency_2d ! MLD-averaged tendency (same units) + + integer :: i,j,k,kp1,n + integer :: tau + real, dimension(isd:ied,jsd:jed) :: mld + real, parameter :: epsln=1.0e-20 ! guard against division by zero at land points + + if (.not.module_is_initialized) then + call mpp_error(FATAL, & + '==>Error from ocean_tracer_diag_mod (compute_budget_mld): module needs initialization') + endif + + tau = Time%tau + + call calc_mixed_layer_depth(Thickness, & + T_prog(index_salt)%field(isd:ied,jsd:jed,:,tau), & + T_prog(index_temp)%field(isd:ied,jsd:jed,:,tau), & + Dens%rho(isd:ied,jsd:jed,:,tau), & + Dens%pressure_at_depth(isd:ied,jsd:jed,:), & + mld(:,:), smooth_mld_input=.false.) + + ! Build the per-layer fractional weight wrk1(i,j,k). + ! Case 1: MLD is shallower than the bottom of the first layer. + wrk1(:,:,:) = 0.0 + k=1 + do j=jsc,jec + do i=isc,iec + if(Grd%tmask(i,j,k)==1.0) then + if(Thickness%depth_zwt(i,j,k) >= mld(i,j)) then + wrk1(i,j,1) = mld(i,j)/Thickness%depth_zwt(i,j,k) + wrk1(i,j,2:nk) = 0.0 + endif + endif + enddo + enddo + + ! Case 2: MLD straddles a layer interface at depth k (k>1). + ! Layers above are fully included; the straddling layer gets a partial weight. + do j=jsc,jec + do i=isc,iec + kloopA: do k=2,nk + if(Grd%tmask(i,j,k)==1.0) then + if(Thickness%depth_zwt(i,j,k) >= mld(i,j) .and. & + Thickness%depth_zwt(i,j,k-1) < mld(i,j)) then + kp1 = min(k+1,nk) + wrk1(i,j,1:k-1) = 1.0 + wrk1(i,j,k) = (mld(i,j)-Thickness%depth_zwt(i,j,k-1))/Thickness%dzt(i,j,k) + wrk1(i,j,kp1:nk) = 0.0 + exit kloopA + endif + endif + enddo kloopA + enddo + enddo + + ! Case 3: MLD is deeper than all model levels — include the full column. + k=nk + do j=jsc,jec + do i=isc,iec + if(Grd%tmask(i,j,k)==1.0) then + if(Thickness%depth_zwt(i,j,k) <= mld(i,j)) then + wrk1(i,j,:) = 1.0 + endif + endif + enddo + enddo + + ! Depth-weighted column sum, then divide by MLD to get depth-average. + wrk1_2d(:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1_2d(i,j) = wrk1_2d(i,j) + wrk1(i,j,k)*tendency(i,j,k) + enddo + enddo + enddo + tendency_2d(:,:) = 0.0 + do j=jsc,jec + do i=isc,iec + if (Grd%tmask(i,j,1)==1.0) then + tendency_2d(i,j) = wrk1_2d(i,j)/(mld(i,j) + epsln) + endif + enddo + enddo + +end subroutine compute_budget_mld +! NAME="compute_budget_mld" + +!####################################################################### +! +! +! +! +! Return the tracer value at the base of the mixed layer by linear +! interpolation between the depth cell centres (depth_zt) that bracket +! the MLD. Three edge cases are handled: +! 1. MLD shallower than the first cell centre: use the surface value. +! 2. MLD between cell centres k and k+1: distance-weighted interpolation +! (W1 = distance above MLD from k, W2 = distance below MLD from k+1; +! value = (tracer_k * W2 + tracer_k+1 * W1) / (W1+W2)). +! 3. MLD deeper than the deepest wet cell: use the bottom value. +! Also returns the MLD itself via mld_out for downstream diagnostics. +! +! Ryan.Holmes, November 2025 +! +! +subroutine compute_tracer_at_mlb(Time, Thickness, Dens, T_prog, tracer, tracer_at_mlb, mld_out) + + type(ocean_time_type), intent(in) :: Time + type(ocean_thickness_type), intent(in) :: Thickness + type(ocean_density_type), intent(in) :: Dens + type(ocean_prog_tracer_type), intent(in) :: T_prog(:) + real, dimension(isd:,jsd:,:), intent(in) :: tracer ! 3D tracer field + real, dimension(isd:,jsd:), intent(inout) :: tracer_at_mlb ! tracer interpolated to MLD base + real, dimension(isd:,jsd:), intent(inout) :: mld_out ! MLD (m), passed out for callers + + integer :: i,j,k,kmt + integer :: tau + real :: W1, W2, denominator_r + real, dimension(isd:ied,jsd:jed) :: mld + real, parameter :: epsln=1.0e-20 ! guard against zero-thickness layers + + if (.not.module_is_initialized) then + call mpp_error(FATAL, & + '==>Error from ocean_tracer_diag_mod (compute_tracer_at_mlb): module needs initialization') + endif + + tau = Time%tau + + call calc_mixed_layer_depth(Thickness, & + T_prog(index_salt)%field(isd:ied,jsd:jed,:,tau), & + T_prog(index_temp)%field(isd:ied,jsd:jed,:,tau), & + Dens%rho(isd:ied,jsd:jed,:,tau), & + Dens%pressure_at_depth(isd:ied,jsd:jed,:), & + mld(:,:), smooth_mld_input=.false.) + + tracer_at_mlb(:,:) = 0.0 + mld_out(:,:) = mld(:,:) + + ! Case 1: MLD shallower than first cell centre — use surface value. + do j=jsc,jec + do i=isc,iec + if(mld(i,j) <= Thickness%depth_zt(i,j,1)) then + tracer_at_mlb(i,j) = tracer(i,j,1) + endif + enddo + enddo + + ! Case 2: MLD between cell centres k and k+1 — linear interpolation. + ! W1 = distance from cell k to MLD; W2 = distance from MLD to cell k+1. + do j=jsc,jec + do i=isc,iec +kloop: do k=1,nk-1 + if(Thickness%depth_zt(i,j,k) < mld(i,j) .and. mld(i,j) <= Thickness%depth_zt(i,j,k+1)) then + if(Grd%tmask(i,j,k+1) > 0) then + W1 = mld(i,j) - Thickness%depth_zt(i,j,k) + W2 = Thickness%depth_zt(i,j,k+1) - mld(i,j) + denominator_r = 1.0/(W1 + W2 + epsln) + tracer_at_mlb(i,j) = (tracer(i,j,k)*W2 + tracer(i,j,k+1)*W1)*denominator_r + exit kloop + endif + endif + enddo kloop + enddo + enddo + + ! Case 3: MLD deeper than the deepest wet cell — use the bottom value. + do j=jsc,jec + do i=isc,iec + kmt = Grd%kmt(i,j) + if(kmt > 0) then + if(mld(i,j) > Thickness%depth_zt(i,j,kmt)) then + tracer_at_mlb(i,j) = tracer(i,j,kmt) + endif + endif + enddo + enddo +end subroutine compute_tracer_at_mlb +! NAME="compute_tracer_at_mlb" !####################################################################### ! diff --git a/src/mom5/ocean_param/lateral/ocean_mixdownslope.F90 b/src/mom5/ocean_param/lateral/ocean_mixdownslope.F90 index eb91042f90..d158699013 100644 --- a/src/mom5/ocean_param/lateral/ocean_mixdownslope.F90 +++ b/src/mom5/ocean_param/lateral/ocean_mixdownslope.F90 @@ -96,6 +96,7 @@ module ocean_mixdownslope_mod use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_sum use ocean_util_mod, only: write_chksum_2d, write_chksum_3d, write_chksum_2d_int use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk1_v implicit none @@ -187,8 +188,9 @@ module ocean_mixdownslope_mod integer :: id_tform_salt_mixdown_on_nrho =-1 -integer, dimension(:), allocatable :: id_mixdownslope -integer, dimension(:), allocatable :: id_mixdownslope_on_nrho +integer, dimension(:), allocatable :: id_mixdownslope +integer, dimension(:), allocatable :: id_mixdownslope_on_nrho +integer, dimension(:), allocatable :: id_mixdownslope_in_mld ! mixdownslope tendency averaged in mixed layer character(len=128) :: version=& @@ -536,8 +538,10 @@ subroutine ocean_mixdownslope_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio allocate (id_mixdownslope(num_prog_tracers)) allocate (id_mixdownslope_on_nrho(num_prog_tracers)) + allocate (id_mixdownslope_in_mld(num_prog_tracers)) id_mixdownslope_on_nrho = -1 id_mixdownslope = -1 + id_mixdownslope_in_mld = -1 do n=1,num_prog_tracers if(T_prog(n)%name == 'temp') then @@ -551,6 +555,11 @@ subroutine ocean_mixdownslope_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio Dens%neutralrho_axes(1:3), Time%model_time, & 'cp*mixdownslope*rho*dzt*temp binned to neutral density', & 'Watt/m^2', missing_value=missing_value, range=(/-1.e20,1.e20/)) + id_mixdownslope_in_mld(n) = register_diag_field ('ocean_model', & + 'mixdownslope_'//trim(T_prog(n)%name)//'_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'cp*mixdownslope*rho*dzt*temp averaged in mixed layer', & + 'Watt/m^3', missing_value=missing_value, range=(/-1.e9,1.e9/)) else id_mixdownslope(n) = register_diag_field ('ocean_model', 'mixdownslope_'//trim(T_prog(n)%name), & Grd%tracer_axes(1:3), Time%model_time, & @@ -561,6 +570,11 @@ subroutine ocean_mixdownslope_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio Dens%neutralrho_axes(1:3), Time%model_time, & 'mixdownslope*rho*dzt*tracer for '//trim(T_prog(n)%name)//' binned to neutral density', & trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e20,1.e20/)) + id_mixdownslope_in_mld(n) = register_diag_field ('ocean_model', & + 'mixdownslope_'//trim(T_prog(n)%name)//'_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'mixdownslope*rho*dzt*tracer averaged in MLD for '//trim(T_prog(n)%name), & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, range=(/-1.e9,1.e9/)) endif enddo @@ -650,6 +664,7 @@ subroutine mixdownslope (Time, Thickness, T_prog, Dens, index_temp, index_salt) real :: mass_sum, tmix real :: mixdownslope_total, mixdownslope_total_r real :: tendency + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld integer :: stdoutunit stdoutunit=stdout() @@ -907,10 +922,10 @@ subroutine mixdownslope (Time, Thickness, T_prog, Dens, index_temp, index_salt) enddo enddo - if(id_mixdownslope(nt) > 0 .or. id_mixdownslope_on_nrho(nt) > 0) then + if(id_mixdownslope(nt) > 0 .or. id_mixdownslope_on_nrho(nt) > 0 .or. id_mixdownslope_in_mld(nt) > 0) then wrk1(:,:,:) = 0.0 do k=1,nk - do j=jsc,jec + do j=jsc,jec do i=isc,iec wrk1(i,j,k) = tend_mix(i,j,k)*mixdownslope_mask(i,j)*T_prog(nt)%conversion enddo @@ -922,6 +937,10 @@ subroutine mixdownslope (Time, Thickness, T_prog, Dens, index_temp, index_salt) if (id_mixdownslope_on_nrho(nt) > 0) then call diagnose_3d_rho(Time, Dens, id_mixdownslope_on_nrho(nt), wrk1) endif + if (id_mixdownslope_in_mld(nt) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_mixdownslope_in_mld(nt), tendency_in_mld(:,:)) + endif endif if(nt==index_temp) then diff --git a/src/mom5/ocean_param/lateral/ocean_sigma_transport.F90 b/src/mom5/ocean_param/lateral/ocean_sigma_transport.F90 index 1f6affe206..ad052ebabf 100644 --- a/src/mom5/ocean_param/lateral/ocean_sigma_transport.F90 +++ b/src/mom5/ocean_param/lateral/ocean_sigma_transport.F90 @@ -208,6 +208,7 @@ module ocean_sigma_transport_mod use ocean_types_mod, only: ocean_options_type, ocean_thickness_type, ocean_density_type use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_sum, write_chksum_2d use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk6 implicit none @@ -237,8 +238,9 @@ module ocean_sigma_transport_mod ! for tracer from sigma diffusion within sigma layer integer, dimension(:), allocatable :: id_sigma_adv_2d ! thickness and density weighted time tendency ! for tracer from sigma advection within sigma layer -integer, dimension(:), allocatable :: id_sigma_smooth ! thickness and density weighted time tendency +integer, dimension(:), allocatable :: id_sigma_smooth ! thickness and density weighted time tendency ! for tracer from smoothing of the sigma thickness +integer, dimension(:), allocatable :: id_sigma_diff_in_mld ! sigma-diffusion tendency averaged in mixed layer integer, dimension(:), allocatable :: id_sigma_diff_xflux ! for i-flux from sigma diffusion integer, dimension(:), allocatable :: id_sigma_diff_yflux ! for j-flux from sigma diffusion integer, dimension(:), allocatable :: id_sigma_diff_xflux_int_z ! for i-flux integrated in z @@ -758,12 +760,14 @@ subroutine ocean_sigma_transport_init(Grid, Domain, Time, Dens, T_prog, Ocean_op allocate (id_sigma_diff_2d(num_prog_tracers)) allocate (id_sigma_adv_2d(num_prog_tracers)) allocate (id_sigma_smooth(num_prog_tracers)) + allocate (id_sigma_diff_in_mld(num_prog_tracers)) id_tracer_sigma = -1 id_sigma_diff = -1 id_sigma_adv = -1 id_sigma_diff_2d = -1 id_sigma_adv_2d = -1 id_sigma_smooth = -1 + id_sigma_diff_in_mld = -1 do n=1,num_prog_tracers @@ -793,7 +797,11 @@ subroutine ocean_sigma_transport_init(Grid, Domain, Time, Dens, T_prog, Ocean_op Grd%tracer_axes(1:2), Time%model_time, & 'thk wghtd sigma-smoothing heating in sigma layer', 'Watts/m^2', & missing_value=missing_value, range=(/-1.e16,1.e16/)) - else + id_sigma_diff_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_sigma_diff_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'thk wghtd sigma-diffusion heating averaged in mixed layer', 'Watts/m^3', & + missing_value=missing_value, range=(/-1.e16,1.e16/)) + else id_sigma_diff(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_sigma_diff', & Grd%tracer_axes(1:3), Time%model_time, & 'thk wghtd sigma-diffusion on '//trim(T_prog(n)%name), trim(T_prog(n)%flux_units),& @@ -815,7 +823,11 @@ subroutine ocean_sigma_transport_init(Grid, Domain, Time, Dens, T_prog, Ocean_op Grd%tracer_axes(1:2), Time%model_time, & 'thk wghtd sigma-smoothing in sigma layer on '//trim(T_prog(n)%name), & trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e16,1.e16/)) - endif + id_sigma_diff_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_sigma_diff_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'thk wghtd sigma-diffusion in MLD on '//trim(T_prog(n)%name), & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, range=(/-1.e16,1.e16/)) + endif enddo @@ -995,6 +1007,7 @@ subroutine sigma_transport (Time, Thickness, Dens, T_prog, Adv_vel, bott_blthick real, dimension(isd:ied,jsd:jed) :: thickness_sigma_old real, dimension(isd:ied,jsd:jed) :: dzt_sigma_tendency real, dimension(isd:ied,jsd:jed) :: tmp + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld real, dimension(2) :: density_test real :: depth_sigma @@ -1468,9 +1481,15 @@ subroutine sigma_transport (Time, Thickness, Dens, T_prog, Adv_vel, bott_blthick call diagnose_2d(Time, Grd, id_tracer_sigma(n), tracer_sigma(:,:,n)) ! send tendency to diagnostic manager - if (id_sigma_diff(n) > 0) then + if (id_sigma_diff(n) > 0) then call diagnose_3d(Time, Grd, id_sigma_diff(n), T_prog(n)%conversion*wrk2(:,:,:)) endif + if (id_sigma_diff_in_mld(n) > 0) then + wrk5(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, wrk2(:,:,:), tendency_in_mld(:,:)) + wrk1(:,:,:) = wrk5(:,:,:) + call diagnose_2d(Time, Grd, id_sigma_diff_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if (id_sigma_adv(n) > 0) then call diagnose_3d(Time, Grd, id_sigma_adv(n), T_prog(n)%conversion*wrk3(:,:,:)) endif diff --git a/src/mom5/ocean_param/lateral/ocean_submesoscale.F90 b/src/mom5/ocean_param/lateral/ocean_submesoscale.F90 index 90aad1e613..0b776d8f04 100644 --- a/src/mom5/ocean_param/lateral/ocean_submesoscale.F90 +++ b/src/mom5/ocean_param/lateral/ocean_submesoscale.F90 @@ -238,6 +238,7 @@ module ocean_submesoscale_mod use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_sum, write_chksum_3d use ocean_tracer_util_mod,only: diagnose_3d_rho use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk1_2d, wrk2_2d, wrk1_v +use ocean_tracer_diag_mod,only: compute_budget_mld implicit none @@ -339,6 +340,7 @@ module ocean_submesoscale_mod integer, dimension(:), allocatable :: id_xflux_submeso_int_z ! vertically integrated i-flux integer, dimension(:), allocatable :: id_yflux_submeso_int_z ! vertically integrated j-flux integer, dimension(:), allocatable :: id_submeso ! tendency from streamfunction portion of submesoscale param +integer, dimension(:), allocatable :: id_submeso_in_mld ! tendency from streamfunction portion of submesoscale param averaged in mixed layer integer, dimension(:), allocatable :: id_submeso_on_nrho ! tendency from streamfunction portion of submesoscale param binned to neutral density integer, dimension(:), allocatable :: id_xflux_subdiff ! i-directed horz diffusive flux @@ -1034,6 +1036,7 @@ subroutine ocean_submesoscale_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio allocate (id_xflux_submeso_int_z(num_prog_tracers)) allocate (id_yflux_submeso_int_z(num_prog_tracers)) allocate (id_submeso(num_prog_tracers)) + allocate (id_submeso_in_mld(num_prog_tracers)) allocate (id_submeso_on_nrho(num_prog_tracers)) allocate (id_xflux_subdiff(num_prog_tracers)) @@ -1077,6 +1080,12 @@ subroutine ocean_submesoscale_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio 'rho*dzt*cp*submesoscale tendency (heating)', & trim(T_prog(n)%flux_units), missing_value=missing_value, & range=(/-1.e10,1.e10/)) + id_submeso_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_submeso_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'rho*dzt*cp*submesoscale tendency (heating) averaged in mixed layer', & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, & + range=(/-1.e10,1.e10/)) id_submeso_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_submeso_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & @@ -1149,6 +1158,12 @@ subroutine ocean_submesoscale_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio 'rho*dzt*submesoscale tendency for '//trim(T_prog(n)%name), & trim(T_prog(n)%flux_units), missing_value=missing_value, & range=(/-1.e10,1.e10/)) + id_submeso_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_submeso_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'rho*dzt*submesoscale tendency averaged in mixed layer for '//trim(T_prog(n)%name), & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, & + range=(/-1.e10,1.e10/)) id_submeso_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_submeso_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & @@ -1265,7 +1280,7 @@ subroutine submeso_restrat(Time, Thickness, Dens, T_prog, surf_blthick) call compute_submeso_skewsion(Thickness, Dens, Time, T_prog) elseif(submeso_advect_flux) then if(submeso_advect_upwind) then - call compute_submeso_upwind(Time, Dens, T_prog) + call compute_submeso_upwind(Thickness, Time, Dens, T_prog) elseif(submeso_advect_sweby) then call compute_submeso_sweby(Thickness, Time, Dens, T_prog) endif @@ -2368,6 +2383,7 @@ subroutine compute_submeso_skewsion(Thickness, Dens, Time, T_prog) type(ocean_prog_tracer_type), intent(inout) :: T_prog(:) real, dimension(isd:ied,jsd:jed) :: eta_tend + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld real :: eta_tend_glob integer :: i,j,k,n,tau @@ -2400,6 +2416,10 @@ subroutine compute_submeso_skewsion(Thickness, Dens, Time, T_prog) if(id_submeso(n) > 0) then call diagnose_3d(Time, Grd, id_submeso(n), T_prog(n)%wrk1(:,:,:)*T_prog(n)%conversion) endif + if (id_submeso_in_mld(n) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, T_prog(n)%wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_submeso_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if(id_submeso_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_submeso_on_nrho(n), T_prog(n)%wrk1*T_prog(n)%conversion) endif @@ -2852,8 +2872,9 @@ end subroutine compute_flux_z ! ! ! -subroutine compute_submeso_upwind(Time, Dens, T_prog) +subroutine compute_submeso_upwind(Thickness, Time, Dens, T_prog) + type(ocean_thickness_type), intent(in) :: Thickness type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens type(ocean_prog_tracer_type), intent(inout) :: T_prog(:) @@ -2862,6 +2883,7 @@ subroutine compute_submeso_upwind(Time, Dens, T_prog) real,dimension(isc:iec,jsc:jec) :: ft2 real,dimension(isc:iec,jsc:jec) :: fe real,dimension(isc:iec,jsc:jec) :: fn + real,dimension(isc:iec,jsc:jec) :: tendency_in_mld integer :: i, j, k, n, kp1, kp2 integer :: tau, taum1 @@ -2985,6 +3007,10 @@ subroutine compute_submeso_upwind(Time, Dens, T_prog) if(id_submeso(n) > 0) then call diagnose_3d(Time, Grd, id_submeso(n), -advect_tendency(:,:,:)*T_prog(n)%conversion) endif + if (id_submeso_in_mld(n) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, -advect_tendency(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_submeso_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if(id_submeso_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_submeso_on_nrho(n), -advect_tendency*T_prog(n)%conversion) endif @@ -3020,6 +3046,7 @@ subroutine compute_submeso_sweby(Thickness, Time, Dens, T_prog) real,dimension(isc:iec,jsc:jec) :: ftp real,dimension(isc:iec,jsc:jec) :: fbt real,dimension(isc:iec,jsc:jec) :: wkm1 + real,dimension(isc:iec,jsc:jec) :: tendency_in_mld integer :: i, j, k, n integer :: kp1, kp2, km1 @@ -3305,7 +3332,10 @@ subroutine compute_submeso_sweby(Thickness, Time, Dens, T_prog) if(id_submeso(n) > 0) then call diagnose_3d(Time, Grd, id_submeso(n), T_prog(n)%wrk1(:,:,:)*T_prog(n)%conversion) endif - + if (id_submeso_in_mld(n) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, T_prog(n)%wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_submeso_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if(id_submeso_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_submeso_on_nrho(n), T_prog(n)%wrk1*T_prog(n)%conversion) endif diff --git a/src/mom5/ocean_param/neutral/ocean_nphysicsC.F90 b/src/mom5/ocean_param/neutral/ocean_nphysicsC.F90 index d0adf3c22e..2b9e06ed2d 100644 --- a/src/mom5/ocean_param/neutral/ocean_nphysicsC.F90 +++ b/src/mom5/ocean_param/neutral/ocean_nphysicsC.F90 @@ -293,6 +293,7 @@ module ocean_nphysicsC_mod use ocean_util_mod, only: write_note, write_line, write_warning use ocean_util_mod, only: diagnose_2d, diagnose_3d use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld use ocean_workspace_mod, only: wrk1_2d, wrk1_v2d, wrk2_v2d use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4 use ocean_workspace_mod, only: wrk1_v, wrk2_v, wrk3_v, wrk4_v @@ -381,6 +382,8 @@ module ocean_nphysicsC_mod integer, dimension(:), allocatable :: id_k33_implicit ! K33 handled implicitly in time integer, dimension(:), allocatable :: id_neutral_physics_ndiffuse_on_nrho ! tendency from neutral diffusion integer, dimension(:), allocatable :: id_neutral_physics_gm_on_nrho ! tendency from GM +integer, dimension(:), allocatable :: id_neutral_physics_ndiffuse_in_mld ! tendency from neutral diffusion +integer, dimension(:), allocatable :: id_neutral_physics_gm_in_mld ! tendency from GM integer, dimension(:), allocatable :: id_flux_x_ndiffuse ! i-directed tracer flux from neutral diffuse integer, dimension(:), allocatable :: id_flux_y_ndiffuse ! j-directed tracer flux from neutral diffuse integer, dimension(:), allocatable :: id_flux_z_ndiffuse ! k-directed tracer flux from neutral diffuse @@ -1137,11 +1140,15 @@ subroutine ocean_nphysicsC_init(Grid, Domain, Time, Time_steps, Thickness, Dens, allocate (id_neutral_physics_gm(num_prog_tracers)) allocate (id_neutral_physics_ndiffuse_on_nrho(num_prog_tracers)) allocate (id_neutral_physics_gm_on_nrho(num_prog_tracers)) + allocate (id_neutral_physics_ndiffuse_in_mld(num_prog_tracers)) + allocate (id_neutral_physics_gm_in_mld(num_prog_tracers)) id_k33_implicit = -1 id_neutral_physics_ndiffuse = -1 id_neutral_physics_gm = -1 id_neutral_physics_ndiffuse_on_nrho = -1 id_neutral_physics_gm_on_nrho = -1 + id_neutral_physics_ndiffuse_in_mld = -1 + id_neutral_physics_gm_in_mld = -1 do n=1,num_prog_tracers id_k33_implicit(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_k33_implicit', & Grd%tracer_axes_wt(1:3), Time%model_time, & @@ -1161,6 +1168,12 @@ subroutine ocean_nphysicsC_init(Grid, Domain, Time, Time_steps, Thickness, Dens, 'rho*dzt*cp*explicit neutral diffusion tendency (heating) binned to neutral density',& trim(T_prog(n)%flux_units), missing_value=missing_value, & range=(/-1.e20,1.e20/)) + id_neutral_physics_ndiffuse_in_mld(n) = register_diag_field ('ocean_model', & + 'neutral_diffusion_in_mld_'//trim(T_prog(n)%name), & + Grd%tracer_axes(1:2), Time%model_time, & + 'rho*dzt*cp*explicit neutral diffusion tendency averaged in mixed layer (heating)',& + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, & + range=(/-1.e10,1.e10/)) id_neutral_physics_gm(n) = register_diag_field ('ocean_model', & 'neutral_gm_'//trim(T_prog(n)%name), & Grd%tracer_axes(1:3), Time%model_time, & @@ -1173,6 +1186,12 @@ subroutine ocean_nphysicsC_init(Grid, Domain, Time, Time_steps, Thickness, Dens, 'rho*dzt*cp*GM stirring (heating) binned to neutral density',& trim(T_prog(n)%flux_units), missing_value=missing_value, & range=(/-1.e20,1.e20/)) + id_neutral_physics_gm_in_mld(n) = register_diag_field ('ocean_model', & + 'neutral_gm_in_mld_'//trim(T_prog(n)%name), & + Grd%tracer_axes(1:2), Time%model_time, & + 'rho*dzt*cp*GM stirring averaged in mixed layer (heating)', & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, & + range=(/-1.e10,1.e10/)) else id_neutral_physics_ndiffuse(n) = register_diag_field ('ocean_model', & 'neutral_diffusion_'//trim(T_prog(n)%name), & @@ -1186,6 +1205,12 @@ subroutine ocean_nphysicsC_init(Grid, Domain, Time, Time_steps, Thickness, Dens, 'rho*dzt*explicit neutral diffusion tendency binned to neutral density for '//trim(T_prog(n)%name),& trim(T_prog(n)%flux_units), missing_value=missing_value, & range=(/-1.e20,1.e20/)) + id_neutral_physics_ndiffuse_in_mld(n) = register_diag_field ('ocean_model', & + 'neutral_diffusion_in_mld_'//trim(T_prog(n)%name), & + Grd%tracer_axes(1:2), Time%model_time, & + 'rho*dzt*explicit neutral diffusion tendency averaged in mixed layer for '//trim(T_prog(n)%name),& + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, & + range=(/-1.e10,1.e10/)) id_neutral_physics_gm(n) = register_diag_field ('ocean_model', & 'neutral_gm_'//trim(T_prog(n)%name), & Grd%tracer_axes(1:3), Time%model_time, & @@ -1198,6 +1223,12 @@ subroutine ocean_nphysicsC_init(Grid, Domain, Time, Time_steps, Thickness, Dens, 'rho*dzt*GM stirring tendency binned to neutral density for '//trim(T_prog(n)%name),& trim(T_prog(n)%flux_units), missing_value=missing_value, & range=(/-1.e20,1.e20/)) + id_neutral_physics_gm_in_mld(n) = register_diag_field ('ocean_model', & + 'neutral_gm_in_mld_'//trim(T_prog(n)%name), & + Grd%tracer_axes(1:2), Time%model_time, & + 'rho*dzt*GM stirring tendency averaged in mixed layer for '//trim(T_prog(n)%name),& + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, & + range=(/-1.e10,1.e10/)) endif enddo @@ -1825,6 +1856,7 @@ subroutine compute_ndiffusion(Time, Thickness, Dens, T_prog) real, dimension(isd:ied,jsd:jed) :: tmp_flux real, dimension(isd:ied,jsd:jed) :: eta_tend + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld real :: eta_tend_glob integer :: i,j,k,n,tau @@ -1884,6 +1916,10 @@ subroutine compute_ndiffusion(Time, Thickness, Dens, T_prog) if(id_neutral_physics_ndiffuse(n) > 0) then call diagnose_3d(Time, Grd, id_neutral_physics_ndiffuse(n), T_prog(n)%wrk1(:,:,:)*T_prog(n)%conversion) endif + if (id_neutral_physics_ndiffuse_in_mld(n) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, T_prog(n)%wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_neutral_physics_ndiffuse_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if(id_neutral_physics_ndiffuse_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_neutral_physics_ndiffuse_on_nrho(n), T_prog(n)%wrk1*T_prog(n)%conversion) endif @@ -1971,6 +2007,7 @@ subroutine compute_gmskewsion(Time, Thickness, Dens, T_prog) real, dimension(isd:ied,jsd:jed) :: tmp_flux real, dimension(isd:ied,jsd:jed) :: eta_tend + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld real :: eta_tend_glob real :: upsilonx, upsilony real :: d_upsilonx_dz, d_upsilony_dz @@ -2033,6 +2070,10 @@ subroutine compute_gmskewsion(Time, Thickness, Dens, T_prog) if(id_neutral_physics_gm(n) > 0) then call diagnose_3d(Time, Grd, id_neutral_physics_gm(n), T_prog(n)%wrk1(:,:,:)*T_prog(n)%conversion) endif + if (id_neutral_physics_gm_in_mld(n) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, T_prog(n)%wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_neutral_physics_gm_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if(id_neutral_physics_gm_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_neutral_physics_gm_on_nrho(n), T_prog(n)%wrk1*T_prog(n)%conversion) endif diff --git a/src/mom5/ocean_param/sources/ocean_rivermix.F90 b/src/mom5/ocean_param/sources/ocean_rivermix.F90 index 1ab7f3848a..9ade62544d 100644 --- a/src/mom5/ocean_param/sources/ocean_rivermix.F90 +++ b/src/mom5/ocean_param/sources/ocean_rivermix.F90 @@ -120,6 +120,7 @@ module ocean_rivermix_mod use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk1_2d use ocean_util_mod, only: diagnose_2d, diagnose_3d, diagnose_sum use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld implicit none @@ -168,6 +169,7 @@ module ocean_rivermix_mod ! for diagnostics integer, dimension(:), allocatable :: id_rivermix +integer, dimension(:), allocatable :: id_rivermix_in_mld integer, dimension(:), allocatable :: id_rivermix_on_nrho integer, dimension(:), allocatable :: id_runoffmix integer, dimension(:), allocatable :: id_calvingmix @@ -538,10 +540,12 @@ subroutine ocean_rivermix_init(Grid, Domain, Time, Time_steps, Dens, & ! register for diag_manager allocate (id_rivermix(num_prog_tracers)) + allocate (id_rivermix_in_mld(num_prog_tracers)) allocate (id_rivermix_on_nrho(num_prog_tracers)) allocate (id_runoffmix(num_prog_tracers)) allocate (id_calvingmix(num_prog_tracers)) id_rivermix = -1 + id_rivermix_in_mld = -1 id_rivermix_on_nrho = -1 id_runoffmix = -1 id_calvingmix = -1 @@ -552,6 +556,10 @@ subroutine ocean_rivermix_init(Grid, Domain, Time, Time_steps, Dens, & Grd%tracer_axes(1:3), Time%model_time, & 'cp*rivermix*rho_dzt*temp', 'Watt/m^2', & missing_value=missing_value, range=(/-1.e10,1.e10/)) + id_rivermix_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_rivermix_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'cp*rivermix*rho_dzt*temp averaged in mixed layer', 'Watt/m^3', & + missing_value=missing_value, range=(/-1.e10,1.e10/)) id_rivermix_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_rivermix_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'cp*rivermix*rho_dzt*temp binned to neutral density', 'Watt/m^2', & @@ -570,6 +578,11 @@ subroutine ocean_rivermix_init(Grid, Domain, Time, Time_steps, Dens, & 'rivermix*rho_dzt*tracer for '//trim(T_prog(n)%name), & trim(T_prog(n)%flux_units), & missing_value=missing_value, range=(/-1.e10,1.e10/)) + id_rivermix_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_rivermix_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'rivermix*rho_dzt*tracer averaged in mixed layer for '//trim(T_prog(n)%name), & + trim(T_prog(n)%flux_units)//'/m', & + missing_value=missing_value, range=(/-1.e10,1.e10/)) id_rivermix_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_rivermix_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'rivermix*rho_dzt*tracer for '//trim(T_prog(n)%name)//' binned to neutral density',& @@ -646,6 +659,7 @@ subroutine rivermix (Time, Thickness, Dens, T_prog, river, runoff, calving, & real, dimension(isd:,jsd:), intent(in) :: calving real, dimension(isd:,jsd:,:,:), intent(inout) :: diff_cbt + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld integer :: i,j,n,tau logical :: river_discharge =.false. logical :: runoff_discharge =.false. @@ -681,6 +695,11 @@ subroutine rivermix (Time, Thickness, Dens, T_prog, river, runoff, calving, & if(id_rivermix(n) > 0) then call diagnose_3d(Time, Grd, id_rivermix(n), T_prog(n)%wrk1(:,:,:)*T_prog(n)%conversion) endif + if (id_rivermix_in_mld(n) > 0) then + tendency_in_mld(:,:) = 0.0 + call compute_budget_mld(Time, Thickness, Dens, T_prog, T_prog(n)%wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_rivermix_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if(id_rivermix_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_rivermix_on_nrho(n), T_prog(n)%wrk1*T_prog(n)%conversion) endif diff --git a/src/mom5/ocean_param/sources/ocean_shortwave.F90 b/src/mom5/ocean_param/sources/ocean_shortwave.F90 index 5cf141da43..e99a3d8349 100644 --- a/src/mom5/ocean_param/sources/ocean_shortwave.F90 +++ b/src/mom5/ocean_param/sources/ocean_shortwave.F90 @@ -47,6 +47,7 @@ module ocean_shortwave_mod use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5 use ocean_util_mod, only: diagnose_2d, diagnose_3d, diagnose_sum use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld implicit none @@ -76,6 +77,7 @@ module ocean_shortwave_mod ! for diagnostics integer :: id_sw_frac =-1 integer :: id_sw_heat =-1 +integer :: id_sw_heat_in_mld =-1 integer :: id_sw_heat_on_nrho=-1 integer :: id_irradiance =-1 @@ -214,6 +216,11 @@ subroutine ocean_shortwave_init(Grid, Domain, Time, Dens, T_prog, vert_coordinat 'W/m^2', missing_value=-1e10, range=(/-1.e10,1.e10/), & standard_name='downwelling_shortwave_flux_in_sea_water') + id_sw_heat_in_mld = register_diag_field ('ocean_model', 'sw_heat_in_mld', & + Grid%tracer_axes(1:2), Time%model_time, 'penetrative shortwave heating averaged in mixed layer', & + 'W/m^3', missing_value=-1e10, range=(/-1.e10,1.e10/), & + standard_name='downwelling_shortwave_flux_in_sea_water') + id_irradiance = register_diag_field ('ocean_model', 'irradiance', & Grid%tracer_axes(1:3), Time%model_time, 'irradiance', & 'W/m^2',missing_value=missing_value, range=(/-1e6,1e6/)) @@ -272,6 +279,7 @@ subroutine sw_source (Time, Thickness, Dens, T_prog, T_diag, swflx, swflx_vis, T real, dimension(isd:,jsd:,:), intent(inout) :: sw_frac_zt real, dimension(isd:,jsd:,:), intent(inout) :: opacity + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld integer :: i,j,k,tau if (.not. use_this_module) return @@ -344,6 +352,10 @@ subroutine sw_source (Time, Thickness, Dens, T_prog, T_diag, swflx, swflx_vis, T call diagnose_3d(Time, Grd, id_sw_frac, sw_frac_zt(:,:,:)) call diagnose_3d(Time, Grd, id_sw_heat, Temp%wrk1(:,:,:)) + if (id_sw_heat_in_mld > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, Temp%wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_sw_heat_in_mld, tendency_in_mld(:,:)) + endif call diagnose_3d(Time, Grd, id_irradiance, T_diag(index_irr)%field(:,:,:)) call watermass_diag(Time, Temp, Dens) diff --git a/src/mom5/ocean_param/vertical/ocean_vert_kpp_mom4p1.F90 b/src/mom5/ocean_param/vertical/ocean_vert_kpp_mom4p1.F90 index d83cb956e5..22c111a865 100644 --- a/src/mom5/ocean_param/vertical/ocean_vert_kpp_mom4p1.F90 +++ b/src/mom5/ocean_param/vertical/ocean_vert_kpp_mom4p1.F90 @@ -263,6 +263,7 @@ module ocean_vert_kpp_mom4p1_mod use ocean_workspace_mod, only: wrk1_2d, wrk2_2d use ocean_util_mod, only: diagnose_2d, diagnose_3d, diagnose_sum use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld implicit none @@ -437,6 +438,7 @@ module ocean_vert_kpp_mom4p1_mod ! for diagnostics integer, dimension(:), allocatable :: id_nonlocal(:) +integer, dimension(:), allocatable :: id_nonlocal_in_mld(:) integer, dimension(:), allocatable :: id_nonlocal_on_nrho(:) integer, dimension(:), allocatable :: id_ghats(:) integer, dimension(:), allocatable :: id_wsfc(:) @@ -873,11 +875,13 @@ subroutine ocean_vert_kpp_mom4p1_init (Grid, Domain, Time, Time_steps, Dens, T_p ! register diagnostics allocate(id_nonlocal(num_prog_tracers)) + allocate(id_nonlocal_in_mld(num_prog_tracers)) allocate(id_nonlocal_on_nrho(num_prog_tracers)) allocate(id_ghats(2)) allocate(id_wsfc(num_prog_tracers)) allocate(id_wbot(num_prog_tracers)) id_nonlocal=-1 + id_nonlocal_in_mld=-1 id_nonlocal_on_nrho=-1 do n = 1, num_prog_tracers if(n==index_temp) then @@ -885,6 +889,10 @@ subroutine ocean_vert_kpp_mom4p1_init (Grid, Domain, Time, Time_steps, Dens, T_p Grd%tracer_axes(1:3), Time%model_time, & 'cp*rho*dzt*nonlocal tendency from KPP', trim(T_prog(n)%flux_units), & missing_value=missing_value, range=(/-1.e10,1.e10/)) + id_nonlocal_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'cp*rho*dzt*nonlocal tendency from KPP averaged in mixed layer', trim(T_prog(n)%flux_units)//'/m', & + missing_value=missing_value, range=(/-1.e10,1.e10/)) id_nonlocal_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'cp*rho*dzt*nonlocal tendency from KPP binned to neutral density', trim(T_prog(n)%flux_units), & @@ -898,6 +906,10 @@ subroutine ocean_vert_kpp_mom4p1_init (Grid, Domain, Time, Time_steps, Dens, T_p Grd%tracer_axes(1:3), Time%model_time, & 'rho*dzt*nonlocal tendency from KPP', trim(T_prog(n)%flux_units), & missing_value=missing_value, range=(/-1.e10,1.e10/)) + id_nonlocal_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, & + 'rho*dzt*nonlocal tendency from KPP averaged in mixed layer', trim(T_prog(n)%flux_units)//'/m', & + missing_value=missing_value, range=(/-1.e10,1.e10/)) id_nonlocal_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_nonlocal_KPP_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'rho*dzt*nonlocal tendency from KPP binned to neutral density', trim(T_prog(n)%flux_units), & @@ -1017,6 +1029,7 @@ subroutine vert_mix_kpp_mom4p1 (aidif, Time, Thickness, Velocity, T_prog, T_diag real, dimension(isd:ied,jsd:jed,nk) :: dbloc1, dbsfc1 real, dimension(isd:ied,jsd:jed) :: frazil + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld real :: smftu, smftv, active_cells real :: temporary, temporary1 integer :: i, j, k, n, ki, ni, kp @@ -1554,6 +1567,10 @@ subroutine vert_mix_kpp_mom4p1 (aidif, Time, Thickness, Velocity, T_prog, T_diag if (id_nonlocal(n) > 0) then call diagnose_3d(Time, Grd, id_nonlocal(n),T_prog(n)%conversion*T_prog(n)%wrk1(:,:,:)) endif + if (id_nonlocal_in_mld(n) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, T_prog(n)%wrk1(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_nonlocal_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if (id_nonlocal_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_nonlocal_on_nrho(n),T_prog(n)%conversion*T_prog(n)%wrk1) endif diff --git a/src/mom5/ocean_param/vertical/ocean_vert_mix.F90 b/src/mom5/ocean_param/vertical/ocean_vert_mix.F90 index fb08116b24..e162f67726 100644 --- a/src/mom5/ocean_param/vertical/ocean_vert_mix.F90 +++ b/src/mom5/ocean_param/vertical/ocean_vert_mix.F90 @@ -276,6 +276,7 @@ module ocean_vert_mix_mod use ocean_util_mod, only: invtri, invtri_bmf, write_timestamp use ocean_util_mod, only: write_chksum_3d, diagnose_2d, diagnose_3d, diagnose_3d_u, diagnose_2d_u, diagnose_sum, diagnose_2d_en use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld use ocean_vert_const_mod, only: ocean_vert_const_init, vert_mix_const use ocean_vert_chen_mod, only: ocean_vert_chen_init, vert_mix_chen, ocean_vert_chen_end use ocean_vert_chen_mod, only: ocean_vert_chen_restart @@ -696,10 +697,14 @@ module ocean_vert_mix_mod integer, dimension(:), allocatable :: id_vdiffuse_impl integer, dimension(:), allocatable :: id_vdiffuse_k33 integer, dimension(:), allocatable :: id_vdiffuse_k33_on_nrho +integer, dimension(:), allocatable :: id_vdiffuse_k33_in_mld integer, dimension(:), allocatable :: id_vdiffuse_diff_cbt +integer, dimension(:), allocatable :: id_vdiffuse_diff_cbt_in_mld integer, dimension(:), allocatable :: id_vdiffuse_diff_cbt_on_nrho integer, dimension(:), allocatable :: id_vdiffuse_diff_cbt_conv +integer, dimension(:), allocatable :: id_vdiffuse_diff_cbt_conv_in_mld integer, dimension(:), allocatable :: id_vdiffuse_sbc +integer, dimension(:), allocatable :: id_vdiffuse_sbc_in_mld integer, dimension(:), allocatable :: id_vdiffuse_sbc_on_nrho integer, dimension(:), allocatable :: id_vdiffuse_bbc integer, dimension(:), allocatable :: id_vdiffuse_diss @@ -991,11 +996,15 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, allocate( id_vdiffuse(num_prog_tracers) ) allocate( id_vdiffuse_impl(num_prog_tracers) ) allocate( id_vdiffuse_k33(num_prog_tracers) ) + allocate( id_vdiffuse_k33_in_mld(num_prog_tracers) ) allocate( id_vdiffuse_k33_on_nrho(num_prog_tracers) ) allocate( id_vdiffuse_diff_cbt(num_prog_tracers) ) + allocate( id_vdiffuse_diff_cbt_in_mld(num_prog_tracers) ) allocate( id_vdiffuse_diff_cbt_on_nrho(num_prog_tracers) ) allocate( id_vdiffuse_diff_cbt_conv(num_prog_tracers) ) + allocate( id_vdiffuse_diff_cbt_conv_in_mld(num_prog_tracers) ) allocate( id_vdiffuse_sbc(num_prog_tracers) ) + allocate( id_vdiffuse_sbc_in_mld(num_prog_tracers) ) allocate( id_vdiffuse_sbc_on_nrho(num_prog_tracers) ) allocate( id_vdiffuse_bbc(num_prog_tracers) ) allocate( id_vdiffuse_diss(num_prog_tracers) ) @@ -1003,11 +1012,15 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, id_vdiffuse =-1 id_vdiffuse_impl =-1 id_vdiffuse_k33 =-1 + id_vdiffuse_k33_in_mld=-1 id_vdiffuse_k33_on_nrho =-1 id_vdiffuse_diff_cbt=-1 + id_vdiffuse_diff_cbt_in_mld=-1 id_vdiffuse_diff_cbt_on_nrho=-1 id_vdiffuse_diff_cbt_conv=-1 + id_vdiffuse_diff_cbt_conv_in_mld=-1 id_vdiffuse_sbc =-1 + id_vdiffuse_sbc_in_mld=-1 id_vdiffuse_sbc_on_nrho=-1 id_vdiffuse_bbc =-1 id_vdiffuse_diss =-1 @@ -1031,6 +1044,10 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, trim(T_prog(n)%name)//'_vdiffuse_k33', Grd%tracer_axes(1:3), & Time%model_time, 'vert diffusion of heat due to K33 from neutral diffusion', 'Watts/m^2',& missing_value=missing_value, range=(/-1.e16,1.e16/)) + id_vdiffuse_k33_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_k33_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'vert diffusion of heat due to K33 from neutral diffusion averaged in mixed layer', 'Watts/m^3',& + missing_value=missing_value, range=(/-1.e16,1.e16/)) id_vdiffuse_k33_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_k33_on_nrho', Dens%neutralrho_axes(1:3), & Time%model_time, 'vert diffusion of heat due to K33 from neutral diffusion binned to neutral density', 'Watts/m^2',& @@ -1039,10 +1056,18 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, trim(T_prog(n)%name)//'_vdiffuse_diff_cbt', Grd%tracer_axes(1:3), & Time%model_time, 'vert diffusion of heat due to diff_cbt', 'Watts/m^2',& missing_value=missing_value, range=(/-1.e16,1.e16/)) + id_vdiffuse_diff_cbt_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'vert diffusion of heat due to diff_cbt averaged in mixed layer', 'Watts/m^3',& + missing_value=missing_value, range=(/-1.e16,1.e16/)) id_vdiffuse_diff_cbt_conv(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_conv', Grd%tracer_axes(1:3), & Time%model_time, 'vert diffusion of heat due to diff_cbt_conv','Watts/m^2', & missing_value=missing_value, range=(/-1.e16,1.e16/)) + id_vdiffuse_diff_cbt_conv_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_conv_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'vert diffusion of heat due to diff_cbt_conv averaged in mixed layer','Watts/m^3', & + missing_value=missing_value, range=(/-1.e16,1.e16/)) id_vdiffuse_diff_cbt_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_on_nrho', Dens%neutralrho_axes(1:3), & Time%model_time, 'vert diffusion of heat due to diff_cbt binned to neutral density', 'Watts/m^2',& @@ -1051,6 +1076,10 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, trim(T_prog(n)%name)//'_vdiffuse_sbc', Grd%tracer_axes(1:3), & Time%model_time, 'vert diffusion of heat due to surface flux', 'Watts/m^2',& missing_value=missing_value, range=(/-1.e16,1.e16/)) + id_vdiffuse_sbc_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_sbc_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'vert diffusion of heat due to surface flux averaged in mixed layer', 'Watts/m^3',& + missing_value=missing_value, range=(/-1.e16,1.e16/)) id_vdiffuse_sbc_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_sbc_on_nrho', Dens%neutralrho_axes(1:3), & Time%model_time, 'vert diffusion of heat due to surface flux binned to neutral density', 'Watts/m^2',& @@ -1080,6 +1109,10 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, trim(T_prog(n)%name)//'_vdiffuse_k33', Grd%tracer_axes(1:3), Time%model_time, & 'vert diffusion due to K33 from neutral diffusion for '//trim(T_prog(n)%longname), & trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1e20,1e20/)) + id_vdiffuse_k33_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_k33_in_mld', Grd%tracer_axes(1:2), Time%model_time, & + 'vert diffusion due to K33 from neutral diffusion averaged in mixed layer for '//trim(T_prog(n)%longname), & + trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1e20,1e20/)) id_vdiffuse_k33_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_k33_on_nrho', Dens%neutralrho_axes(1:3), & Time%model_time, 'vert diffusion due to K33 from neutral diffusion '//trim(T_prog(n)%longname)//' binned to neutral density',& @@ -1088,6 +1121,10 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, trim(T_prog(n)%name)//'_vdiffuse_diff_cbt', Grd%tracer_axes(1:3), Time%model_time, & 'vert diffusion due to diff_cbt for '//trim(T_prog(n)%longname), & trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1e20,1e20/)) + id_vdiffuse_diff_cbt_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_in_mld', Grd%tracer_axes(1:2), Time%model_time, & + 'vert diffusion due to diff_cbt averaged in mixed layer for '//trim(T_prog(n)%longname), & + trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1e20,1e20/)) id_vdiffuse_diff_cbt_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time, & 'vert diffusion due to diff_cbt for '//trim(T_prog(n)%longname)//' binned to neutral density',& @@ -1096,10 +1133,18 @@ subroutine ocean_vert_mix_init (Grid, Domain, Time, Dens, Velocity, Time_steps, trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_conv', Grd%tracer_axes(1:3),Time%model_time, & 'vert diffusion due to diff_cbt_conv for '//trim(T_prog(n)%longname), & trim(T_prog(n)%flux_units), missing_value=missing_value,range=(/-1e20,1e20/)) + id_vdiffuse_diff_cbt_conv_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_diff_cbt_conv_in_mld', Grd%tracer_axes(1:2),Time%model_time, & + 'vert diffusion due to diff_cbt_conv averaged in mixed layer for '//trim(T_prog(n)%longname), & + trim(T_prog(n)%flux_units), missing_value=missing_value,range=(/-1e20,1e20/)) id_vdiffuse_sbc(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_sbc', Grd%tracer_axes(1:3), Time%model_time,& 'vert diffusion due to surface flux for '//trim(T_prog(n)%longname), & trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1e20,1e20/)) + id_vdiffuse_sbc_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_vdiffuse_sbc_in_mld', Grd%tracer_axes(1:2), Time%model_time,& + 'vert diffusion due to surface flux averaged in mixed layer for '//trim(T_prog(n)%longname), & + trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1e20,1e20/)) id_vdiffuse_sbc_on_nrho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_vdiffuse_sbc_on_nrho', Dens%neutralrho_axes(1:3), Time%model_time,& 'vert diffusion due to surface flux for '//trim(T_prog(n)%longname)//' binned to neutral density',& @@ -5234,8 +5279,9 @@ subroutine vert_diffuse_implicit_diag(Time, Thickness, Dens, T_prog, diff_cbt, w integer :: taup1, tau integer :: i, j, k, kp, kp1, nmix - real :: dTdz, diffusivity - + real :: dTdz, diffusivity + real,dimension(isd:ied,jsd:jed) :: tendency_in_mld + real,dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d tau = Time%tau taup1 = Time%taup1 @@ -5278,7 +5324,7 @@ subroutine vert_diffuse_implicit_diag(Time, Thickness, Dens, T_prog, diff_cbt, w ! impacts from surface boundary fluxes - if(id_vdiffuse_sbc(n) > 0 .or. id_vdiffuse_sbc_on_nrho(n) > 0) then + if(id_vdiffuse_sbc(n) > 0 .or. id_vdiffuse_sbc_in_mld(n) > 0 .or. id_vdiffuse_sbc_on_nrho(n) > 0) then wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 wrk1(:,:,:) = 0.0 @@ -5292,6 +5338,12 @@ subroutine vert_diffuse_implicit_diag(Time, Thickness, Dens, T_prog, diff_cbt, w if (id_vdiffuse_sbc(n) > 0) then call diagnose_3d(Time, Grd, id_vdiffuse_sbc(n), wrk1(:,:,:)*T_prog(n)%conversion) endif + if (id_vdiffuse_sbc_in_mld(n) > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_vdiffuse_sbc_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if (id_vdiffuse_sbc_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_vdiffuse_sbc_on_nrho(n), wrk1*T_prog(n)%conversion) endif @@ -5314,7 +5366,7 @@ subroutine vert_diffuse_implicit_diag(Time, Thickness, Dens, T_prog, diff_cbt, w ! impacts from diff_cbt - if(id_vdiffuse_diff_cbt(n) > 0 .or. id_vdiffuse_diff_cbt_on_nrho(n) > 0) then + if(id_vdiffuse_diff_cbt(n) > 0 .or. id_vdiffuse_diff_cbt_in_mld(n) > 0 .or. id_vdiffuse_diff_cbt_on_nrho(n) > 0) then wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 wrk1(:,:,:) = 0.0 @@ -5334,13 +5386,19 @@ subroutine vert_diffuse_implicit_diag(Time, Thickness, Dens, T_prog, diff_cbt, w if (id_vdiffuse_diff_cbt(n) > 0) then call diagnose_3d(Time, Grd, id_vdiffuse_diff_cbt(n), wrk1(:,:,:)*T_prog(n)%conversion) endif + if (id_vdiffuse_diff_cbt_in_mld(n) > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_vdiffuse_diff_cbt_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if (id_vdiffuse_diff_cbt_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_vdiffuse_diff_cbt_on_nrho(n), wrk1*T_prog(n)%conversion) endif endif ! impacts from diff_cbt_conv - if(id_vdiffuse_diff_cbt_conv(n) > 0) then + if(id_vdiffuse_diff_cbt_conv(n) > 0 .or. id_vdiffuse_diff_cbt_in_mld(n) > 0) then wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 wrk1(:,:,:) = 0.0 @@ -5357,11 +5415,19 @@ subroutine vert_diffuse_implicit_diag(Time, Thickness, Dens, T_prog, diff_cbt, w enddo enddo enddo - call diagnose_3d(Time, Grd, id_vdiffuse_diff_cbt_conv(n),wrk1(:,:,:)*T_prog(n)%conversion) + if (id_vdiffuse_diff_cbt_conv(n) > 0) then + call diagnose_3d(Time, Grd, id_vdiffuse_diff_cbt_conv(n),wrk1(:,:,:)*T_prog(n)%conversion) + endif + if (id_vdiffuse_diff_cbt_conv_in_mld(n) > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_vdiffuse_diff_cbt_conv_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif endif ! impacts from K33 - if(id_vdiffuse_k33(n) > 0 .or. id_vdiffuse_k33_on_nrho(n) > 0) then + if(id_vdiffuse_k33(n) > 0 .or. id_vdiffuse_k33_in_mld(n) > 0 .or. id_vdiffuse_k33_on_nrho(n) > 0) then wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 wrk1(:,:,:) = 0.0 @@ -5381,6 +5447,12 @@ subroutine vert_diffuse_implicit_diag(Time, Thickness, Dens, T_prog, diff_cbt, w if (id_vdiffuse_k33(n) > 0) then call diagnose_3d(Time, Grd, id_vdiffuse_k33(n), wrk1(:,:,:)*T_prog(n)%conversion) endif + if (id_vdiffuse_k33_in_mld(n) > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_vdiffuse_k33_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if (id_vdiffuse_k33_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_vdiffuse_k33_on_nrho(n), wrk1*T_prog(n)%conversion) endif diff --git a/src/mom5/ocean_tracers/ocean_frazil.F90 b/src/mom5/ocean_tracers/ocean_frazil.F90 index 6a248b9e8e..09ff0cd115 100644 --- a/src/mom5/ocean_tracers/ocean_frazil.F90 +++ b/src/mom5/ocean_tracers/ocean_frazil.F90 @@ -127,6 +127,7 @@ module ocean_frazil_mod use ocean_types_mod, only: ocean_prog_tracer_type, ocean_diag_tracer_type, ocean_density_type use ocean_workspace_mod, only: wrk1_2d, wrk1 use ocean_util_mod, only: diagnose_2d, diagnose_3d, diagnose_sum +use ocean_tracer_diag_mod,only: compute_budget_mld #if defined(ACCESS_CM) || defined(ACCESS_OM) use auscom_ice_parameters_mod, only: pop_icediag, do_ice @@ -166,6 +167,7 @@ module ocean_frazil_mod logical :: used integer :: id_frazil_2d=-1 integer :: id_frazil_3d=-1 +integer :: id_frazil_3d_in_mld=-1 integer :: id_frazil_3d_int_z=-1 integer :: id_total_ocean_frazil=-1 integer :: id_temp_freeze=-1 @@ -481,6 +483,9 @@ subroutine ocean_frazil_init (Domain, Grid, Time, Time_steps, Ocean_options, & id_frazil_3d = register_diag_field ('ocean_model', 'frazil_3d', Grd%tracer_axes(1:3), & Time%model_time, 'ocn frazil heat flux over time step', 'W/m^2', & missing_value=missing_value, range=(/-1.e10,1.e10/)) + id_frazil_3d_in_mld = register_diag_field ('ocean_model', 'frazil_3d_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'ocn frazil heat flux over time step averaged in mixed layer', 'W/m^3', & + missing_value=missing_value, range=(/-1.e10,1.e10/)) id_frazil_3d_int_z = register_diag_field ('ocean_model', 'frazil_3d_int_z', Grd%tracer_axes(1:2), & Time%model_time, 'Vertical sum of ocn frazil heat flux over time step', 'W/m^2', & missing_value=missing_value, range=(/-1.e10,1.e10/)) @@ -531,6 +536,8 @@ subroutine compute_frazil_heating (Time, Thickness, Dens, T_prog, T_diag) real :: tf_num, tf_den, tfreeze real :: s, sqrts real :: press + real,dimension(isd:ied,jsd:jed) :: tendency_in_mld + real,dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d if(.not. use_this_module) return @@ -557,6 +564,14 @@ subroutine compute_frazil_heating (Time, Thickness, Dens, T_prog, T_diag) Time%model_time, rmask=Grd%tmask(:,:,:), & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk) endif + if (id_frazil_3d_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = T_diag(index_frazil)%field(:,:,:)*dtimer + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + used = send_data(id_frazil_3d_in_mld, tendency_in_mld(:,:), & + Time%model_time, rmask=Grd%tmask(:,:,1), & + is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) + endif if (id_frazil_3d_int_z > 0 .or. id_total_ocean_frazil > 0) then wrk1_2d(:,:) = 0.0 @@ -724,6 +739,12 @@ subroutine compute_frazil_heating (Time, Thickness, Dens, T_prog, T_diag) if (id_frazil_3d > 0) then call diagnose_3d(Time, Grd, id_frazil_3d, T_diag(index_frazil)%field(:,:,:)*dtimer) endif + if (id_frazil_3d_in_mld > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = T_diag(index_frazil)%field(:,:,:)*dtimer + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_frazil_3d_in_mld, tendency_in_mld(:,:)) + endif if (id_frazil_3d_int_z > 0 .or. id_total_ocean_frazil > 0) then wrk1_2d(:,:) = 0.0 do k=1,nk diff --git a/src/mom5/ocean_tracers/ocean_tracer.F90 b/src/mom5/ocean_tracers/ocean_tracer.F90 index 208c91752e..d580484c20 100644 --- a/src/mom5/ocean_tracers/ocean_tracer.F90 +++ b/src/mom5/ocean_tracers/ocean_tracer.F90 @@ -205,7 +205,7 @@ module ocean_tracer_mod use ocean_tpm_mod, only: ocean_tpm_init use ocean_tpm_util_mod, only: otpm_set_tracer_package use ocean_tracer_advect_mod, only: horz_advect_tracer, vert_advect_tracer -use ocean_tracer_diag_mod, only: send_tracer_variance +use ocean_tracer_diag_mod, only: send_tracer_variance, compute_budget_mld, compute_tracer_at_mlb use ocean_tracer_util_mod, only: rebin_onto_rho, diagnose_mass_of_layer use ocean_tracer_util_mod, only: tracer_prog_chksum, tracer_diag_chksum, tracer_min_max use ocean_thickness_mod, only: update_E_thickness @@ -314,6 +314,7 @@ module ocean_tracer_mod logical :: used integer, allocatable, dimension(:) :: id_eta_smooth integer, allocatable, dimension(:) :: id_eta_smooth_on_nrho +integer, allocatable, dimension(:) :: id_eta_smooth_in_mld integer, allocatable, dimension(:) :: id_pbot_smooth integer, allocatable, dimension(:) :: id_prog integer, allocatable, dimension(:) :: id_progT @@ -327,6 +328,9 @@ module ocean_tracer_mod integer, allocatable, dimension(:) :: id_tendency_concL integer, allocatable, dimension(:) :: id_tendency_concT integer, allocatable, dimension(:) :: id_tendency +integer, allocatable, dimension(:) :: id_tendency_in_mld +integer, allocatable, dimension(:) :: id_tracer_in_mld +integer, allocatable, dimension(:) :: id_tracer_at_mlb integer, allocatable, dimension(:) :: id_tendency_on_nrho integer, allocatable, dimension(:) :: id_tendencyL integer, allocatable, dimension(:) :: id_tendencyT @@ -766,6 +770,7 @@ function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, T allocate( T_prog (num_prog_tracers) ) allocate( id_eta_smooth (num_prog_tracers) ) allocate( id_eta_smooth_on_nrho(num_prog_tracers) ) + allocate( id_eta_smooth_in_mld(num_prog_tracers) ) allocate( id_pbot_smooth (num_prog_tracers) ) allocate( id_prog (num_prog_tracers) ) allocate( id_prog_explicit (num_prog_tracers) ) @@ -774,6 +779,9 @@ function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, T allocate( id_prog_on_depth (num_prog_tracers) ) allocate( id_tendency_conc (num_prog_tracers) ) allocate( id_tendency (num_prog_tracers) ) + allocate( id_tendency_in_mld(num_prog_tracers) ) + allocate( id_tracer_in_mld(num_prog_tracers) ) + allocate( id_tracer_at_mlb(num_prog_tracers) ) allocate( id_tendency_on_nrho(num_prog_tracers) ) allocate( id_tendency_expl (num_prog_tracers) ) allocate( id_surf_tracer (num_prog_tracers) ) @@ -801,6 +809,7 @@ function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, T id_eta_smooth(:) = -1 id_eta_smooth_on_nrho(:)= -1 + id_eta_smooth_in_mld(:)= -1 id_pbot_smooth(:) = -1 id_prog(:) = -1 id_prog_explicit(:) = -1 @@ -809,6 +818,9 @@ function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, T id_prog_on_depth(:) = -1 id_tendency_conc(:) = -1 id_tendency(:) = -1 + id_tendency_in_mld(:)= -1 + id_tracer_in_mld(:)= -1 + id_tracer_at_mlb(:)= -1 id_tendency_on_nrho(:)= -1 id_tendency_expl(:) = -1 id_surf_tracer(:) = -1 @@ -1401,11 +1413,28 @@ function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, T trim(prog_name)//'_tendency', Grd%tracer_axes(1:3), & Time%model_time, 'time tendency for tracer '//trim(prog_longname), & trim(T_prog(n)%flux_units), missing_value=missing_value, range=range_array) + id_tendency_in_mld(n) = register_diag_field ('ocean_model', & + trim(prog_name)//'_tendency_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'time tendency averged in mixed layer for tracer '//trim(prog_longname), & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, range=range_array) + id_tracer_in_mld(n) = register_diag_field ('ocean_model', & + trim(prog_name)//'_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'tracer averaged in mixed layer * rho for tracer '//trim(prog_longname), & + trim(T_prog(n)%units)//' kg m-3', missing_value=missing_value, range=range_array) + id_tracer_at_mlb(n) = register_diag_field ('ocean_model', & + trim(prog_name)//'_at_mlb', Grd%tracer_axes(1:2), & + Time%model_time, 'tracer at base of mixed layer for tracer '//trim(prog_longname), & + trim(T_prog(n)%units), missing_value=missing_value, range=range_array) id_eta_smooth(n) = register_diag_field ('ocean_model', & trim(prog_name)//'_eta_smooth', Grd%tracer_axes(1:2), & Time%model_time, 'surface smoother for ' // trim(prog_name), & trim(T_prog(n)%flux_units), & missing_value=missing_value, range=range_array) + id_eta_smooth_in_mld(n) = register_diag_field ('ocean_model', & + trim(prog_name)//'_eta_smooth_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'surface smoother averaged in mixed layer for ' // trim(prog_name), & + trim(T_prog(n)%flux_units)//'/m', & + missing_value=missing_value, range=range_array) id_pbot_smooth(n) = register_diag_field ('ocean_model', & trim(prog_name)//'_pbot_smooth', Grd%tracer_axes(1 :2), & Time%model_time, 'bottom smoother for ' // trim(prog_name), & @@ -1428,11 +1457,28 @@ function ocean_prog_tracer_init (Grid, Thickness, Ocean_options, Domain, Time, T trim(prog_name)//'_tendency', Grd%tracer_axes(1:3), & Time%model_time, 'time tendency for tracer '//trim(prog_longname), & trim(T_prog(n)%flux_units), missing_value=missing_value) + id_tendency_in_mld(n) = register_diag_field ('ocean_model', & + trim(prog_name)//'_tendency_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'time tendency averaged in mixed layer for tracer '//trim(prog_longname), & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value) + id_tracer_in_mld(n) = register_diag_field ('ocean_model', & + trim(prog_name)//'_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'tracer averaged in mixed layer * rho for tracer '//trim(prog_longname), & + trim(T_prog(n)%units)//' kg m-3', missing_value=missing_value) + id_tracer_at_mlb(n) = register_diag_field ('ocean_model', & + trim(prog_name)//'_at_mlb', Grd%tracer_axes(1:2), & + Time%model_time, 'tracer at base of mixed layer for tracer '//trim(prog_longname), & + trim(T_prog(n)%units), missing_value=missing_value) id_eta_smooth(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_eta_smooth', Grd%tracer_axes(1:2), & Time%model_time, 'surface smoother for ' // trim(T_prog(n)%name), & trim(T_prog(n)%flux_units), & missing_value=missing_value) + id_eta_smooth_in_mld(n) = register_diag_field ('ocean_model', & + trim(T_prog(n)%name)//'_eta_smooth_in_mld', Grd%tracer_axes(1:2), & + Time%model_time, 'surface smoother averaged in mixed layer for ' // trim(T_prog(n)%name), & + trim(T_prog(n)%flux_units)//'/m', & + missing_value=missing_value) id_pbot_smooth(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_pbot_smooth', Grd%tracer_axes(1:2), & Time%model_time, 'bottom smoother for ' // trim(T_prog(n)%name), & @@ -2266,6 +2312,9 @@ subroutine update_ocean_tracer (Time, Dens, Adv_vel, Thickness, pme, diff_cbt, & integer :: i, j, k, kbot, n integer :: taum1, tau, taup1 + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d + integer :: stdoutunit stdoutunit=stdout() @@ -2394,6 +2443,13 @@ subroutine update_ocean_tracer (Time, Dens, Adv_vel, Thickness, pme, diff_cbt, & if (id_eta_smooth(n)> 0) then call diagnose_2d(Time, Grd, id_eta_smooth(n), T_prog(n)%eta_smooth(:,:)*T_prog(n)%conversion) endif + if (id_eta_smooth_in_mld(n)> 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = 0.0 + tendency_3d(:,:,1) = T_prog(n)%eta_smooth(:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_eta_smooth_in_mld(n), tendency_in_mld(:,:)*T_prog(n)%conversion) + endif if (id_eta_smooth_on_nrho(n)> 0) then wrk1(:,:,:) = 0.0 k=1 @@ -3853,6 +3909,8 @@ subroutine send_tracer_diagnostics(Time, T_prog, T_diag, Thickness, Dens, use_bl type(ocean_density_type), intent(in) :: Dens logical, intent(in) :: use_blobs + real, dimension(isd:ied,jsd:jed) :: tendency_in_mld + real, dimension(isd:ied,jsd:jed,1:nk) :: tendency_3d integer :: i,j,k,kbot,n integer :: taum1,tau,taup1 real :: total_tracer @@ -3942,7 +4000,7 @@ subroutine send_tracer_diagnostics(Time, T_prog, T_diag, Thickness, Dens, use_bl endif ! time tendency for tracer mass per horizontal area - if (id_tendency(n) > 0 .or. id_tendency_on_nrho(n) > 0) then + if (id_tendency(n) > 0 .or. id_tendency_in_mld(n) > 0 .or. id_tendency_on_nrho(n) > 0) then wrk1(:,:,:) = 0.0 do k=1,nk do j=jsc,jec @@ -3957,11 +4015,39 @@ subroutine send_tracer_diagnostics(Time, T_prog, T_diag, Thickness, Dens, use_bl if (id_tendency(n) > 0) then call diagnose_3d(Time, Grd, id_tendency(n),wrk1(:,:,:)) endif + if (id_tendency_in_mld(n) > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_tendency_in_mld(n), tendency_in_mld(:,:)) + endif if (id_tendency_on_nrho(n) > 0) then call diagnose_3d_rho(Time, Dens, id_tendency_on_nrho(n),wrk1) endif endif + if (id_tracer_in_mld(n) > 0) then + wrk1(:,:,:) = 0.0 + do k=1,nk + do j=jsc,jec + do i=isc,iec + wrk1(i,j,k) = T_prog(n)%field(i,j,k,tau)*Thickness%rho_dzt(i,j,k,tau) + enddo + enddo + enddo + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = wrk1(:,:,:) + call compute_budget_mld(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_tracer_in_mld(n), tendency_in_mld(:,:)) + endif + + if (id_tracer_at_mlb(n) > 0) then + tendency_in_mld(:,:) = 0.0 + tendency_3d(:,:,:) = T_prog(n)%field(:,:,:,tau) + call compute_tracer_at_mlb(Time, Thickness, Dens, T_prog, tendency_3d(:,:,:), tendency_in_mld(:,:), wrk1_2d(:,:)) + call diagnose_2d(Time, Grd, id_tracer_at_mlb(n), tendency_in_mld(:,:)) + endif + ! time tendency for tracer concentration do k=1,nk do j=jsd,jed diff --git a/src/mom5/ocean_tracers/ocean_tracer_advect.F90 b/src/mom5/ocean_tracers/ocean_tracer_advect.F90 index 116ef68810..419baa0420 100644 --- a/src/mom5/ocean_tracers/ocean_tracer_advect.F90 +++ b/src/mom5/ocean_tracers/ocean_tracer_advect.F90 @@ -170,6 +170,7 @@ module ocean_tracer_advect_mod use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk1_2d use ocean_util_mod, only: diagnose_2d, diagnose_3d use ocean_tracer_util_mod, only: diagnose_3d_rho +use ocean_tracer_diag_mod, only: compute_budget_mld implicit none @@ -330,8 +331,8 @@ module ocean_tracer_advect_mod integer :: index_temp_sq=-1 integer :: index_salt_sq=-1 -character(len=256) :: version='CVS $Id: ocean_tracer_advect.F90,v 20.0 2013/12/14 00:17:22 fms Exp $' -character(len=256) :: tagname='Tag $Name: tikal $' +character(len=256) :: version='CVS $Id: ocean_tracer_advect.F90,v 20.0 2013/12/14 00:17:22 fms Exp $' +character(len=256) :: tagname='Tag $Name: tikal $' type(ocean_domain_type), pointer :: Dom =>NULL() @@ -374,6 +375,7 @@ module ocean_tracer_advect_mod integer :: id_tform_salt_advect_on_nrho integer, dimension(:), allocatable :: id_tracer_advection +integer, dimension(:), allocatable :: id_tracer_advection_in_mld integer, dimension(:), allocatable :: id_tracer_advection_on_nrho integer, dimension(:), allocatable :: id_tracer2_advection integer, dimension(:), allocatable :: id_tracer_adv_diss @@ -775,6 +777,7 @@ subroutine advection_diag_init (Time, Dens, T_prog) stdoutunit=stdout() allocate (id_tracer_advection(num_prog_tracers)) + allocate (id_tracer_advection_in_mld(num_prog_tracers)) allocate (id_tracer_advection_on_nrho(num_prog_tracers)) allocate (id_tracer2_advection(num_prog_tracers)) allocate (id_tracer_adv_diss(num_prog_tracers)) @@ -792,6 +795,7 @@ subroutine advection_diag_init (Time, Dens, T_prog) allocate (id_yflux_adv_int_z(num_prog_tracers)) id_tracer_advection =-1 + id_tracer_advection_in_mld =-1 id_tracer_advection_on_nrho =-1 id_tracer2_advection=-1 id_tracer_adv_diss =-1 @@ -814,6 +818,9 @@ subroutine advection_diag_init (Time, Dens, T_prog) id_tracer_advection(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_advection', & Grd%tracer_axes(1:3), Time%model_time, 'cp*rho*dzt*advection tendency', & trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e18,1.e18/)) + id_tracer_advection_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_advection_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, 'cp*rho*dzt/mld*advection tendency averaged in mixed layer', & + trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, range=(/-1.e18,1.e18/)) id_tracer_advection_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_advection_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, 'cp*rho*dzt*advection tendency binned to neutral density',& trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e18,1.e18/)) @@ -872,6 +879,9 @@ subroutine advection_diag_init (Time, Dens, T_prog) id_tracer_advection(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_advection', & Grd%tracer_axes(1:3), Time%model_time, 'rho*dzt*advection tendency', & 'kg/(sec*m^2)', missing_value=missing_value, range=(/-1.e18,1.e18/)) + id_tracer_advection_in_mld(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_advection_in_mld', & + Grd%tracer_axes(1:2), Time%model_time, 'rho*dzt/mld*advection tendency averaged in mixed layer', & + 'kg/(sec*m^3)', missing_value=missing_value, range=(/-1.e18,1.e18/)) id_tracer_advection_on_nrho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_advection_on_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, 'rho*dzt*advection tendency binned to neutral density',& 'kg/(sec*m^2)', missing_value=missing_value, range=(/-1.e18,1.e18/)) @@ -2105,6 +2115,7 @@ subroutine vert_advect_tracer(Time, Adv_vel, Dens, Thickness, T_prog, Tracer, nt integer :: tau, taum1 integer :: i,j,k + real,dimension(isd:ied,jsd:jed) :: tendency_in_mld if(zero_tracer_advect_vert) return @@ -2179,6 +2190,10 @@ subroutine vert_advect_tracer(Time, Adv_vel, Dens, Thickness, T_prog, Tracer, nt if(id_tracer_advection(ntracer) > 0) then call diagnose_3d(Time, Grd, id_tracer_advection(ntracer), Tracer%conversion*advect_tendency(:,:,:)) endif + if(id_tracer_advection_in_mld(ntracer) > 0) then + call compute_budget_mld(Time, Thickness, Dens, T_prog, advect_tendency(:,:,:), tendency_in_mld(:,:)) + call diagnose_2d(Time, Grd, id_tracer_advection_in_mld(ntracer), Tracer%conversion*tendency_in_mld(:,:)) + endif if(id_tracer_advection_on_nrho(ntracer) > 0) then call diagnose_3d_rho(Time, Dens, id_tracer_advection_on_nrho(ntracer), Tracer%conversion*advect_tendency) endif