module ocean_tracer_diag_mod ! ! S.M. Griffies ! ! ! ! Routines for tracer diagnostics ! ! ! ! Routines for tracer diagnostics. Some are printed to ascii output, some are sent ! to diagnostic manager. ! ! ! ! ! Number of days between which compute the tracer conservation diagnostics. ! ! ! Number of time steps between which compute the diagnostics. ! ! ! ! Set true for help with debugging the diagnostic for mixing. ! ! ! Set true for more help with debugging the diagnostic for mixing. ! Lots of output. ! ! ! Set true for more help with debugging the diagnostic for mixing. ! Lots of output. ! ! ! Set true for more help with debugging the diagnostic for mixing. ! Lots of output. ! ! ! Number of 1-2-1 smooths applied to kappa_sort ! ! ! Number of 1-2-1 smooths applied to rho_sort ! ! ! min vertical density gradient (kg/m^3/m) used in computing kappa sorted ! in the diagnostic mixing sorted. ! ! ! max vertical density gradient (kg/m^3/m) used in computing kappa sorted ! ! ! Critical buoyancy difference relative to surface for computing mixed ! layer depth. Default buoyancy_crit=0.0003. ! ! ! Critical temperature difference relative to surface for computing ! mixed_layer_depth_dtheta . Default dtheta_crit=2.0. ! ! ! Days over which time average the thickness weighted density before taking its ! time tendency for use in computing the effective diapycnal diffusivity. ! ! ! ! The preTEOS10 EOS used in MOM requires salinity to ! use the Practical Salinity Scale (pss). This scale is ! also known as the Practical Salinity Unit (psu). ! ! However, salinity as an absolute concentration in ! parts per thousand is more convenient to use when ! performing budget analyses such as in this module. ! Conversion between pss and ppt depends on the precise ! ratio of ions in the seawater. Hence, the conversion ! is not constant. However, it is close to a constant, ! as reported in Jackett etal (2004). For purposes of ! budgets, we take this conversion as a constant. ! The conversion is ! ! s(ppt) = psu2ppt * s(psu) ! ! where again s(psu) is what MOM carries as its ! prognostic salinity field when preTEOS10 EOS is used. ! ! Jackett etal (2004), correcting a type in equation (53) ! of Feistel (2003), report that ! ! s(ppt) = 1.004867 * s(psu) ! ! ! ! ! Smooth the diagnosed mixed layer depth. Default smooth_mld=.false. ! ! ! ! Smooth the diagnosed mixed layer depth to be used for subduction ! diagnostics. Default smooth_mld_for_subduction=.true. ! ! ! ! Set true to do bitwise exact global sum. When it is false, the global ! sum will be non-bitwise_exact, but will significantly increase efficiency. ! The default value is false. ! ! ! ! use constants_mod, only: epsln, c2dbars use diag_manager_mod, only: register_diag_field, send_data, need_data use fms_mod, only: open_namelist_file, check_nml_error, close_file use fms_mod, only: write_version_number, FATAL, WARNING, stdout, stdlog use mpp_domains_mod, only: mpp_global_field, mpp_global_sum, mpp_global_max, mpp_global_min use mpp_domains_mod, only: mpp_update_domains use mpp_domains_mod, only: BITWISE_EXACT_SUM, NON_BITWISE_EXACT_SUM use mpp_mod, only: input_nml_file, mpp_error, mpp_max use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE use time_manager_mod, only: time_type, increment_time use ocean_density_mod, only: density_delta_sfc use ocean_domains_mod, only: get_local_indices use ocean_obc_mod, only: ocean_obc_tracer_flux, ocean_obc_mass_flux use ocean_operators_mod, only: FAX, FAY, BAX, BAY, FMX, FMY, FDX_T, FDY_T, S2D use ocean_parameters_mod, only: DEPTH_BASED, TWO_LEVEL, THREE_LEVEL, missing_value use ocean_parameters_mod, only: rho0, rho0r, sec_in_yr, sec_in_yr_r, grav use ocean_parameters_mod, only: onehalf, onefourth use ocean_tracer_util_mod, only: tracer_min_max, dzt_min_max, sort_shell_array use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_velocity_type use ocean_types_mod, only: ocean_prog_tracer_type, ocean_diag_tracer_type use ocean_types_mod, only: ocean_external_mode_type, ocean_density_type, ocean_adv_vel_type use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type, ocean_thickness_type use ocean_types_mod, only: ocean_lagrangian_type use ocean_util_mod, only: write_timestamp, diagnose_2d use ocean_tracer_util_mod, only: diagnose_3d_rho use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4 use ocean_workspace_mod, only: wrk1_2d, wrk2_2d, wrk3_2d, wrk4_2d use ocean_workspace_mod, only: wrk1_v, wrk2_v , wrk3_v use ocean_workspace_mod, only: wrk1_v2d implicit none private #include real :: dtts real :: dteta real :: dtime real :: dtimer real :: aidif integer :: num_prog_tracers integer :: num_diag_tracers logical :: have_obc logical :: use_blobs integer :: index_temp integer :: index_salt integer :: index_frazil ! for area of domain real :: cellarea real :: cellarea_r ! for vertical coordinate integer :: vert_coordinate_class=1 ! for diagnostics clocks integer :: id_compute_subduction integer :: id_compute_tracer_mld integer :: id_mixed_layer_depth integer :: id_mixed_layer_depth_dtheta integer :: id_potrho_mixed_layer integer :: id_diagnose_depth_of_potrho integer :: id_diagnose_depth_of_theta integer :: id_diagnose_tracer_on_rho integer :: id_diagnose_tracer_zrho_on_rho integer :: id_tracer_numerical integer :: id_diagnose_mixing integer :: id_conservation integer :: id_total_tracer integer :: id_total_mass integer :: id_total_volume integer :: id_tracer_integrals integer :: id_tracer_change integer :: id_tracer_land_cell_check type(ocean_grid_type), pointer :: Grd =>NULL() type(ocean_domain_type), pointer :: Dom =>NULL() logical :: module_is_initialized = .FALSE. integer :: tendency=0 character(len=128) :: version=& '$Id: ocean_tracer_diag.F90,v 20.0 2013/12/14 00:12:55 fms Exp $' character (len=128) :: tagname = & '$Name: tikal $' ! for tracer conservation: need to know num_prog_tracers to allocate real, dimension(:,:,:), allocatable :: tracer_source real, dimension(:,:,:), allocatable :: tracer_runoff real, dimension(:,:,:), allocatable :: tracer_calving real, dimension(:,:,:), allocatable :: tracer_otf real, dimension(:,:,:), allocatable :: tracer_stf real, dimension(:,:,:), allocatable :: tracer_btf real, dimension(:,:,:), allocatable :: tracer_pme real, dimension(:,:), allocatable :: tracer_frazil real, dimension(:,:,:), allocatable :: tracer_tend real, dimension(:,:,:), allocatable :: tracer_eta_smooth real, dimension(:,:,:), allocatable :: tracer_pbot_smooth real, dimension(:,:) , allocatable :: pme_flux real, dimension(:,:) , allocatable :: runoff_flux real, dimension(:,:) , allocatable :: source_flux real, dimension(:,:) , allocatable :: calving_flux real, dimension(:,:) , allocatable :: eta_smooth real, dimension(:,:) , allocatable :: pbot_smooth real, dimension(:,:) , allocatable :: area_t real, dimension(:), allocatable :: tracer_start real, dimension(:), allocatable :: tracer_final real :: mass_start=0.0 real :: mass_final=0.0 real :: volume_start=0.0 real :: volume_final=0.0 real :: tracer_conserve_days=30.0 ! number of days for tracer conservation integer :: itts_tracer, itte_tracer ! itt starting and ending tracer conservation integer :: itts_mass, itte_mass ! itt starting and ending mass conservation ! for saving out the total tracer content and mass integer, allocatable, dimension(:) :: id_prog_total integer, allocatable, dimension(:) :: id_prog_surface_ave integer, allocatable, dimension(:) :: id_prog_surface_area_ave integer, allocatable, dimension(:) :: id_prog_global_ave integer :: id_mass_seawater=-1 integer :: id_volume_seawater=-1 integer :: id_press_ave=-1 ! for diagnosing mixing between temperature classes. logical :: debug_diagnose_mixingA=.false. logical :: debug_diagnose_mixingB=.false. logical :: debug_diagnose_mixingC=.false. logical :: debug_diagnose_mixingD=.false. real :: rho_grad_max=1e28 !max vertical density gradient (kg/m^3/m) used in computing kappa real :: rho_grad_min=1e-5 !min vertical density gradient (kg/m^3/m) used in computing kappa real :: alpha=1.0 !should=alpha_linear_eos, but choose 1.0 to increase precision real :: thick_column(2) !thickness of sorted column real :: mass_column(2) !mass of sorted column integer :: smooth_kappa_sort=0 !for smoothing the diagnosed diffusivity integer :: nsortpts ! number points to be sorted = number wet tracer points integer, dimension(:),allocatable :: nlayer ! # wet points per layer when specify this as constant real, dimension(:), allocatable :: wetarea ! horizontal area as function of depth real, dimension(:), allocatable :: normalize !1.0/(rho0*wetarea(k)) real, dimension(:), allocatable :: deriv_lay ! vertical density deriv at bottom of tracer cell real, dimension(:), allocatable :: flux_lay ! dianeutral tracer flux centered at bottom of sorted tracer cell real, dimension(:), allocatable :: kappa_lay ! diagnosed dianeutral diffusivity (m^2/sec) in layer real, dimension(:,:), allocatable :: rho_lay ! sorted density (kg/m^3) averaged onto layers real, dimension(:,:), allocatable :: dzt_rho_lay ! dzt*rho (kg/m^2) of a layer real, dimension(:,:), allocatable :: dzt_lay ! thickness (m) of a density layer real, dimension(:,:), allocatable :: dzw_lay ! distance (m) between density centers real, dimension(:,:), allocatable :: mass_lay ! mass (kg) of a layer real, dimension(:,:), allocatable :: volume_lay ! volume (m^3) of a layer ! for the spurious mixing diagnostics integer :: id_kappa_simple =-1 integer :: id_kappa_sort =-1 integer :: id_temp_sort =-1 ! frequency which compute certain ascii diagnostics integer :: diag_step = -1 ! for diagnostic manager logical :: used ! for output integer :: unit=6 ! for mixed layer depth integer :: id_mld =-1 integer :: id_mld_sq =-1 integer :: id_mld_nrho =-1 ! for subduction diagnostics integer :: id_subduction =-1 integer :: id_subduction_mld =-1 integer :: id_subduction_dhdt =-1 integer :: id_subduction_horz =-1 integer :: id_subduction_vert =-1 integer :: id_subduction_nrho =-1 integer :: id_subduction_dhdt_nrho =-1 integer :: id_subduction_horz_nrho =-1 integer :: id_subduction_vert_nrho =-1 integer, dimension(:), allocatable :: id_subduction_tracer integer, dimension(:), allocatable :: id_subduction_dhdt_tracer integer, dimension(:), allocatable :: id_subduction_horz_tracer integer, dimension(:), allocatable :: id_subduction_vert_tracer integer, dimension(:), allocatable :: id_subduction_mld_zflux_diff logical :: compute_subduction_diags = .false. ! for tracer integrated over mixed layer integer, dimension(:), allocatable :: id_tracer_mld logical :: compute_tracer_mld_diags = .false. ! for mixed layer depth based solely on depth where SST - temp(k) = dtheta integer :: id_mld_dtheta =-1 ! for depth of isopycnal and potential temperature surfaces integer :: id_depth_of_potrho=-1 integer :: id_depth_of_theta =-1 ! for tracer_on_rho integer, dimension(:), allocatable :: id_tracer_on_rho integer, dimension(:), allocatable :: id_tracer_zrho_on_rho integer, dimension(:), allocatable :: id_global_variance integer, dimension(:), allocatable :: id_k_variance ! for potential density mixed layer integer :: id_potrho_mix_depth=-1 integer :: id_potrho_mix_base=-1 real :: buoyancy_crit=0.0003 ! (m/s^2) critical buoyancy difference relative to surface ! for potential temperature based mixed layer real :: dtheta_crit=2.0 ! (degC) critical temperature difference relative to surface ! for salt budgets real :: psu2ppt=1.004867 ! conversion from psu to ppt concentration ! for bitwise exact global sums independent of PE number logical :: do_bitwise_exact_sum = .false. integer :: global_sum_flag ! for computation of frazil contributions to heat, we need to have ! frazil_factor equal to that used in the nml from ocean_frazil_mod. ! If tendency=TWO_LEVEL and if use GFDL SIS ocean model, then frazil_factor=1.0 ! If tendency=THREE_LEVEL and if use GFDL SIS ocean model, then frazil_factor=.50 real :: frazil_factor=1.0 ! for smoothing the diagnosed mld logical :: smooth_mld=.false. ! for mld used in subduction diagnostic logical :: smooth_mld_for_subduction=.true. public ocean_tracer_diag_init public ocean_tracer_diagnostics public calc_mixed_layer_depth public calc_potrho_mixed_layer public send_tracer_variance public diagnose_eta_tend_3dflux private compute_subduction private compute_tracer_mld private mixed_layer_depth private mixed_layer_depth_dtheta private potrho_mixed_layer private tracer_change private tracer_integrals private tracer_land_cell_check private tracer_conservation private mass_conservation private diagnose_kappa_simple private diagnose_kappa_sort private diagnose_depth_of_potrho private diagnose_depth_of_theta private diagnose_tracer_on_rho private diagnose_tracer_zrho_on_rho private total_tracer private klevel_total_tracer private total_mass private total_volume private klevel_total_mass private send_total_tracer private send_global_ave_tracer private send_global_ave_pressure private send_surface_ave_tracer private send_surface_area_ave_tracer private send_total_mass private send_total_volume namelist /ocean_tracer_diag_nml/ tracer_conserve_days , diag_step, psu2ppt, & debug_diagnose_mixingA, debug_diagnose_mixingB, & debug_diagnose_mixingC, debug_diagnose_mixingD, & smooth_kappa_sort, rho_grad_min, rho_grad_max, & buoyancy_crit, do_bitwise_exact_sum, frazil_factor, & smooth_mld, smooth_mld_for_subduction, dtheta_crit contains !####################################################################### ! ! ! ! Initialize the ocean_tracer_diag module containing subroutines ! diagnosing tracer related properties of the simulation. These are ! not terms in the equations, but rather they are diagnosed from ! terms. ! ! subroutine ocean_tracer_diag_init(Grid, Domain, Time, Time_steps, Thickness, T_prog, T_diag, Dens, & ver_coordinate_class, blobs, obc) type(ocean_grid_type), target, intent(in) :: Grid type(ocean_domain_type), target, intent(in) :: Domain type(ocean_time_type), intent(in) :: Time type(ocean_time_steps_type), intent(in) :: Time_steps type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: T_prog(:) type(ocean_diag_tracer_type), intent(in) :: T_diag(:) type(ocean_density_type), intent(inout) :: Dens integer, intent(in) :: ver_coordinate_class logical, intent(in) :: blobs logical, intent(in) :: obc integer :: n, ioun, io_status, ierr, tau integer :: stdoutunit,stdlogunit stdoutunit=stdout();stdlogunit=stdlog() if (module_is_initialized) return module_is_initialized = .TRUE. num_prog_tracers = size(T_prog,1) num_diag_tracers = size(T_diag,1) call write_version_number(version, tagname) #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=ocean_tracer_diag_nml, iostat=io_status) ierr = check_nml_error(io_status,'ocean_tracer_diag_nml') #else ioun = open_namelist_file() read(ioun, ocean_tracer_diag_nml, iostat=io_status) ierr = check_nml_error(io_status,'ocean_tracer_diag_nml') call close_file(ioun) #endif write (stdlogunit, ocean_tracer_diag_nml) write (stdoutunit,'(/)') write (stdoutunit, ocean_tracer_diag_nml) if (diag_step == 0) diag_step = 1 if(do_bitwise_exact_sum) then global_sum_flag = BITWISE_EXACT_SUM else global_sum_flag = NON_BITWISE_EXACT_SUM endif vert_coordinate_class = ver_coordinate_class #ifndef MOM_STATIC_ARRAYS call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec) ni = Grid%ni nj = Grid%nj nk = Grid%nk #endif Dom => Domain Grd => Grid dtts = Time_steps%dtts dteta = Time_steps%dteta dtime = Time_steps%dtime_t tendency = Time_steps%tendency aidif = Time_steps%aidif dtimer = 1.0/(dtime+epsln) have_obc = obc use_blobs = blobs tau = Time%tau 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 index_frazil=-1 do n=1,num_diag_tracers if(T_diag(n)%name == 'frazil' ) index_frazil = n enddo cellarea = Grd%tcellsurf cellarea_r = 1.0/(epsln + Grd%tcellsurf) allocate(area_t(isd:ied,jsd:jed) ) area_t(:,:) = Grd%dat(:,:)*Grd%tmask(:,:,1) if(have_obc) area_t(:,:) = area_t(:,:)*Grd%obc_tmask(:,:) if(tendency==THREE_LEVEL) then write (stdoutunit,'(/a)') & 'Note: tracer and mass/volume conservation tests based on time_tendency==threelevel.' elseif(tendency==TWO_LEVEL) then write (stdoutunit,'(/a)') & 'Note: tracer and mass/volume conservation tests based on time_tendency==twolevel.' endif write (stdoutunit,'(/a,f5.2,a)') & ' Note: Set frazil_factor = ',frazil_factor,' for computation of heat diagnostics.' write (stdoutunit,'(a)') & ' Be sure this agrees with the value set in nml for ocean_frazil_mod' ! registers to diagnostic manager id_mass_seawater = register_diag_field ('ocean_model', 'total_mass_seawater', Time%model_time,& 'total mass of liquid seawater', 'kg', missing_value=missing_value, & range=(/-1e1,1e25/), standard_name='sea_water_mass') id_volume_seawater = register_diag_field ('ocean_model', 'total_volume_seawater', Time%model_time,& 'total volume of liquid seawater', 'm^3', missing_value=missing_value, & range=(/-1e1,1e25/), standard_name='sea_water_volume') allocate( id_prog_total(num_prog_tracers) ) allocate( id_prog_surface_ave(num_prog_tracers) ) allocate( id_prog_surface_area_ave(num_prog_tracers) ) allocate( id_prog_global_ave(num_prog_tracers) ) id_prog_total = -1 id_prog_surface_ave = -1 id_prog_surface_area_ave= -1 id_prog_global_ave = -1 allocate(id_tracer_on_rho(num_prog_tracers)) allocate(id_tracer_zrho_on_rho(num_prog_tracers)) id_tracer_on_rho=-1 allocate(id_global_variance(num_prog_tracers)) allocate(id_k_variance(num_prog_tracers)) id_global_variance(:) = -1 id_k_variance(:) = -1 allocate( id_subduction_tracer(num_prog_tracers) ) allocate( id_subduction_dhdt_tracer(num_prog_tracers) ) allocate( id_subduction_horz_tracer(num_prog_tracers) ) allocate( id_subduction_vert_tracer(num_prog_tracers) ) allocate( id_subduction_mld_zflux_diff(num_prog_tracers) ) allocate( id_tracer_mld(num_prog_tracers) ) id_subduction_tracer(:) = -1 id_subduction_dhdt_tracer(:) = -1 id_subduction_horz_tracer(:) = -1 id_subduction_vert_tracer(:) = -1 id_subduction_mld_zflux_diff(:) = -1 id_tracer_mld(:) = -1 do n=1,num_prog_tracers if(T_prog(n)%name == 'temp') then id_prog_total(n) = register_diag_field ('ocean_model','total_ocean_heat', & Time%model_time, 'Total heat in the liquid ocean referenced to 0degC','Joule/1e25',& missing_value=missing_value, range=(/0.0,1e20/)) id_tracer_on_rho(n) = register_diag_field ('ocean_model', 'temp_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & 'temperature on potential density surface', 'C', & missing_value=missing_value, range=(/-5.0,100./)) id_tracer_zrho_on_rho(n) = register_diag_field ('ocean_model', 'temp_zrho_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & 'temperature*dz/drho on potential density surface', 'degC * m/(kg/m^3)', & missing_value=missing_value, range=(/-1.e9,1.e9/)) id_prog_global_ave(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_ave',& Time%model_time, 'Global mean '//trim(T_prog(n)%name)//' in liquid seawater' & ,trim(T_prog(n)%units), missing_value=missing_value, range=(/-10.0,1000.0/), & standard_name='sea_water_potential_temperature') id_global_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_variance',& Time%model_time, 'Global '//trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/), & standard_name='global_temperature_variance') id_k_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_variance', & Grd%tracer_axes(3:3), Time%model_time, trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/), & standard_name='temperature_variance') elseif(T_prog(n)%name =='salt') then id_prog_total(n) = register_diag_field ('ocean_model', 'total_ocean_salt', & Time%model_time, 'total mass of salt in liquid seawater', & 'kg/1e18', missing_value=missing_value, range=(/-1e2,1e10/)) id_prog_global_ave(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_ave',& Time%model_time, 'Global mean '//trim(T_prog(n)%name)//' in liquid seawater' & ,trim(T_prog(n)%units), missing_value=missing_value, range=(/-10.0,1000.0/), & standard_name='sea_water_salinity') id_tracer_on_rho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & trim(T_prog(n)%name)//' on potential density surface', trim(T_prog(n)%units), & missing_value=missing_value, range=(/0.0,1.e3/)) id_tracer_zrho_on_rho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_zrho_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & trim(T_prog(n)%name)//' *dz/drho on potential density surface', & trim(T_prog(n)%units)//'*m/(kg/m^3)', & missing_value=missing_value, range=(/-1.e9,1.e9/)) id_global_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_variance',& Time%model_time, 'Global '//trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/), & standard_name='global_salinity_variance') id_k_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_variance', & Grd%tracer_axes(3:3), Time%model_time, trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/), & standard_name='salt_variance') elseif(T_prog(n)%name(1:3) =='age') then id_prog_total(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_total', & Time%model_time, 'mass integrated '//trim(T_prog(n)%name), & 'yr', missing_value=missing_value, range=(/-1e2,1e10/)) id_prog_global_ave(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_ave',& Time%model_time, 'Global mean '//trim(T_prog(n)%name)//' in liquid seawater' & ,trim(T_prog(n)%units), missing_value=missing_value, range=(/-10.0,1000.0/)) id_tracer_on_rho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & trim(T_prog(n)%name)//' on potential density surface', trim(T_prog(n)%units), & missing_value=missing_value, range=(/0.0,1.e3/)) id_tracer_zrho_on_rho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_zrho_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & trim(T_prog(n)%name)//' *dz/drho on potential density surface', & trim(T_prog(n)%units)//'*m/(kg/m^3)', & missing_value=missing_value, range=(/-1.e9,1.e9/)) id_global_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_variance',& Time%model_time, 'Global '//trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/)) id_k_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_variance', & Grd%tracer_axes(3:3), Time%model_time, trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/)) else id_prog_total(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_total', & Time%model_time, 'mass integrated '//trim(T_prog(n)%name), & 'kg/1e18', missing_value=missing_value, range=(/-1e6,1e6/)) id_prog_global_ave(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_ave', & Time%model_time, 'Global mass weighted mean '//trim(T_prog(n)%name)//' in liquid seawater' & ,trim(T_prog(n)%units), missing_value=missing_value, range=(/-10.0,1000.0/)) id_tracer_on_rho(n) = register_diag_field ('ocean_model', trim(T_prog(n)%name)//'_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & trim(T_prog(n)%name)//' on potential density surface', trim(T_prog(n)%units), & missing_value=missing_value, range=(/0.0,1.e3/)) id_tracer_zrho_on_rho(n) = register_diag_field ('ocean_model', & trim(T_prog(n)%name)//'_zrho_on_rho', & Dens%potrho_axes(1:3), Time%model_time, & trim(T_prog(n)%name)//' *dz/drho on potential density surface', & trim(T_prog(n)%units)//'*m/(kg/m^3)', & missing_value=missing_value, range=(/-1.e9,1.e9/)) id_global_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_global_variance',& Time%model_time, 'Global '//trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/)) id_k_variance(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_variance', & Grd%tracer_axes(3:3), Time%model_time, trim(T_prog(n)%name)//' variance', & trim(T_prog(n)%units), missing_value=missing_value, range=(/-1e20,1e20/)) endif id_prog_surface_ave(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_surface_ave', & Time%model_time, 'Global mass weighted mean surface '//trim(T_prog(n)%name)//' in liquid seawater'& ,trim(T_prog(n)%units), missing_value=missing_value, range=(/-10.0,1000.0/)) id_prog_surface_area_ave(n) = register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_surface_area_ave',& Time%model_time, 'Global area weighted mean surface '//trim(T_prog(n)%name)//' in liquid seawater' & ,trim(T_prog(n)%units), missing_value=missing_value, range=(/-10.0,1000.0/)) id_subduction_tracer(n) = & register_diag_field ('ocean_model','subduction_'//trim(T_prog(n)%name), & Grd%tracer_axes(1:2), Time%model_time, & 'Tracer transport across mld from subduction for '//trim(T_prog(n)%name),& 'kg/sec * '//trim(T_prog(n)%units), missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_tracer(n) > 0) compute_subduction_diags = .true. id_subduction_dhdt_tracer(n) = & register_diag_field ('ocean_model','subduction_dhdt_'//trim(T_prog(n)%name),& Grd%tracer_axes(1:2), Time%model_time, & 'Tracer transport across mld from dhdt for '//trim(T_prog(n)%name), & 'kg/sec * '//trim(T_prog(n)%units), missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_dhdt_tracer(n) > 0) compute_subduction_diags = .true. id_subduction_horz_tracer(n) = & register_diag_field ('ocean_model','subduction_horz_'//trim(T_prog(n)%name), & Grd%tracer_axes(1:2),Time%model_time, & 'Tracer transport across mld from from horz advection for '//trim(T_prog(n)%name),& 'kg/sec * '//trim(T_prog(n)%units), missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_horz_tracer(n) > 0) compute_subduction_diags = .true. id_subduction_vert_tracer(n) = & register_diag_field ('ocean_model','subduction_vert_'//trim(T_prog(n)%name), & Grd%tracer_axes(1:2), Time%model_time, & 'Tracer transport across mld from vertical advection for '//trim(T_prog(n)%name),& 'kg/sec * '//trim(T_prog(n)%units), missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_vert_tracer(n) > 0) compute_subduction_diags = .true. id_subduction_mld_zflux_diff(n) = & register_diag_field ('ocean_model','subduction_mld_zflux_diff'//trim(T_prog(n)%name),& Grd%tracer_axes(1:2), Time%model_time, & 'Tracer transport from vertical diffusion at base of mld for '//trim(T_prog(n)%name),& 'kg/sec * '//trim(T_prog(n)%units), missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_mld_zflux_diff(n) > 0) compute_subduction_diags = .true. id_tracer_mld(n) = & register_diag_field ('ocean_model',trim(T_prog(n)%name)//'_mld', & Grd%tracer_axes(1:2), Time%model_time, & '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', & Grd%tracer_axes(3:3), Time%model_time, & 'kappa from sorting', 'm^2/sec', & missing_value=missing_value, range=(/-1e6,1e6/)) id_temp_sort = register_diag_field ('ocean_model','temp_sort', & Grd%tracer_axes(3:3), Time%model_time, & 'sorted temp', 'C', & missing_value=missing_value, range=(/-1e6,1e6/)) id_kappa_simple = register_diag_field ('ocean_model','kappa_simple', & Grd%tracer_axes(3:3), Time%model_time, & 'kappa from horz avg', 'm^2/sec', & missing_value=missing_value, range=(/-1e6,1e6/)) id_mld_dtheta = register_diag_field ('ocean_model', 'mld_dtheta', & Grd%tracer_axes(1:2), Time%model_time, & 'mixed layer depth determined by temperature criteria ', 'm',& missing_value=missing_value, range=(/0.0,1.e6/)) id_mld = register_diag_field ('ocean_model', 'mld', & Grd%tracer_axes(1:2), Time%model_time, & 'mixed layer depth determined by density criteria ', 'm',& missing_value=missing_value, range=(/0.0,1.e6/), & standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t') id_mld_sq = register_diag_field ('ocean_model', 'mld_sq', & Grd%tracer_axes(1:2), Time%model_time, & 'squared mixed layer depth determined by density criteria', 'm^2',& missing_value=missing_value, range=(/0.0,1.e12/), & standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t') id_mld_nrho = register_diag_field ('ocean_model', 'mld_nrho', & Grd%tracer_axes(1:2), Time%model_time, & 'neutral density at base of mixed layer depth', 'kg/m^3',& missing_value=missing_value, range=(/0.0,1.e6/)) id_subduction = register_diag_field ('ocean_model', 'subduction', & Grd%tracer_axes(1:2), Time%model_time, & 'rate of mass transferred below the mixed layer base', & 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction > 0) compute_subduction_diags = .true. id_subduction_mld = register_diag_field ('ocean_model', 'subduction_mld',& Grd%tracer_axes(1:2), Time%model_time, & 'mixed layer depth used for subduction diagnostics', & 'm', missing_value=missing_value, range=(/-1.e2,1.e20/)) if(id_subduction_mld > 0) compute_subduction_diags = .true. id_subduction_dhdt = register_diag_field ('ocean_model', 'subduction_dhdt', & Grd%tracer_axes(1:2), Time%model_time, & 'rate of mass transferred below the mixed layer base due to dh/dt', & 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_dhdt > 0) compute_subduction_diags = .true. id_subduction_horz = register_diag_field ('ocean_model', 'subduction_horz', & Grd%tracer_axes(1:2), Time%model_time, & 'rate of mass transferred below the mixed layer base due to horz advect',& 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_horz > 0) compute_subduction_diags = .true. id_subduction_vert = register_diag_field ('ocean_model', 'subduction_vert', & Grd%tracer_axes(1:2), Time%model_time, & 'rate of mass transferred below the mixed layer base due to vert advect',& 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_vert > 0) compute_subduction_diags = .true. id_subduction_nrho = register_diag_field ('ocean_model', 'subduction_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'rate of mass transferred below the mixed layer base as binned to neutral rho classes',& 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_nrho > 0) compute_subduction_diags = .true. id_subduction_dhdt_nrho = register_diag_field ('ocean_model', 'subduction_dhdt_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'rate of mass transferred below the mixed layer base due to dhdt as binned to neutral rho classes',& 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_dhdt_nrho > 0) compute_subduction_diags = .true. id_subduction_horz_nrho = register_diag_field ('ocean_model', 'subduction_horz_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'rate of mass transferred below the mixed layer base due to horz velocity as binned to neutral rho classes',& 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_horz_nrho > 0) compute_subduction_diags = .true. id_subduction_vert_nrho = register_diag_field ('ocean_model', 'subduction_vert_nrho', & Dens%neutralrho_axes(1:3), Time%model_time, & 'rate of mass transferred below the mixed layer base due to vert velocity as binned to neutral rho classes',& 'kg/sec', missing_value=missing_value, range=(/-1.e20,1.e20/)) if(id_subduction_vert_nrho > 0) compute_subduction_diags = .true. id_depth_of_potrho = register_diag_field ('ocean_model', 'depth_of_potrho', & Dens%potrho_axes(1:3), Time%model_time, & 'depth of potential density surface', 'm', & missing_value=missing_value, range=(/0.0,1.e10/)) id_depth_of_theta = register_diag_field ('ocean_model', 'depth_of_theta', & Dens%theta_axes(1:3), Time%model_time, & 'depth of potential temp surface', 'm', & missing_value=missing_value, range=(/0.0,1.e10/)) id_potrho_mix_depth = register_diag_field ('ocean_model','potrho_mix_depth', & Grd%tracer_axes(1:2), Time%model_time, & 'Depth of potential density mixed layer','m', & missing_value=missing_value, range=(/-1e6,1e6/)) id_potrho_mix_base = register_diag_field ('ocean_model','potrho_mix_base', & Grd%tracer_axes(1:2), Time%model_time, & 'Potential density at mixed layer base','kg/m^3', & missing_value=missing_value, range=(/-1e6,1e6/)) id_press_ave = register_diag_field('ocean_model','press_ave', & Time%model_time, 'global mean absolute ocean pressure, including surface loading',& 'dbar', missing_value=missing_value, range=(/-1e1,1e20/)) ! define clock integers id_compute_subduction = mpp_clock_id('(Ocean tracer_diag: compute_subduction)',grain=CLOCK_ROUTINE) id_compute_tracer_mld = mpp_clock_id('(Ocean tracer_diag: compute_tracer_mld)',grain=CLOCK_ROUTINE) id_mixed_layer_depth_dtheta = mpp_clock_id('(Ocean tracer_diag: mld_dtheta)' ,grain=CLOCK_ROUTINE) id_mixed_layer_depth = mpp_clock_id('(Ocean tracer_diag: mld)' ,grain=CLOCK_ROUTINE) id_potrho_mixed_layer = mpp_clock_id('(Ocean tracer_diag: rho_mld)' ,grain=CLOCK_ROUTINE) id_diagnose_depth_of_potrho = mpp_clock_id('(Ocean tracer_diag: potrho depth)' ,grain=CLOCK_ROUTINE) id_diagnose_depth_of_theta = mpp_clock_id('(Ocean tracer_diag: theta depth)' ,grain=CLOCK_ROUTINE) id_diagnose_tracer_on_rho = mpp_clock_id('(Ocean tracer_diag: tracer on rho)' ,grain=CLOCK_ROUTINE) id_diagnose_tracer_zrho_on_rho = mpp_clock_id('(Ocean tracer_diag: tracer_zrho on rho)',grain=CLOCK_ROUTINE) id_tracer_numerical = mpp_clock_id('(Ocean tracer_diag: numerical)' ,grain=CLOCK_ROUTINE) id_diagnose_mixing = mpp_clock_id('(Ocean tracer_diag: diag mixing' ,grain=CLOCK_ROUTINE) id_conservation = mpp_clock_id('(Ocean tracer_diag: conserve)' ,grain=CLOCK_ROUTINE) id_total_mass = mpp_clock_id('(Ocean tracer_diag: total water mass)' ,grain=CLOCK_ROUTINE) id_total_volume = mpp_clock_id('(Ocean tracer_diag: total water volume)',grain=CLOCK_ROUTINE) id_total_tracer = mpp_clock_id('(Ocean tracer_diag: total tracer)' ,grain=CLOCK_ROUTINE) id_tracer_integrals = mpp_clock_id('(Ocean tracer_diag: integrals)' ,grain=CLOCK_ROUTINE) id_tracer_change = mpp_clock_id('(Ocean tracer_diag: change)' ,grain=CLOCK_ROUTINE) id_tracer_land_cell_check = mpp_clock_id('(Ocean tracer_diag: land check)' ,grain=CLOCK_ROUTINE) ! compute tau value for Dens%mld_subduction 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,:), & Dens%mld_subduction(:,:), smooth_mld_for_subduction) end subroutine ocean_tracer_diag_init ! NAME="ocean_tracer_diag_init" !####################################################################### ! ! ! ! Call diagnostics related to the tracer fields. ! subroutine ocean_tracer_diagnostics(Time, Thickness, T_prog, T_diag, Dens, & Ext_mode, Velocity, Adv_vel, & diff_cbt, pme, melt, runoff, calving) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: T_prog(:) type(ocean_diag_tracer_type), intent(in) :: T_diag(:) type(ocean_density_type), intent(inout) :: Dens type(ocean_external_mode_type), intent(in) :: Ext_mode type(ocean_velocity_type), intent(in) :: Velocity type(ocean_adv_vel_type), intent(in) :: Adv_vel real, dimension(isd:,jsd:,:,:), intent(in) :: diff_cbt real, dimension(isd:,jsd:), intent(in) :: pme real, dimension(isd:,jsd:), intent(in) :: melt real, dimension(isd:,jsd:), intent(in) :: runoff real, dimension(isd:,jsd:), intent(in) :: calving type(time_type) :: next_time integer :: tau, taum1, taup1 integer :: n if(size(T_prog,1) /= num_prog_tracers) then call mpp_error(FATAL,& '==>Error from ocean_tracer_diag_mod (ocean_tracer_diagnostics): T_prog has wrong dimensions') endif tau = Time%tau taum1 = Time%taum1 taup1 = Time%taup1 next_time = increment_time(Time%model_time, int(dtts), 0) ! compute mixed layer depth and send to diagnostics manager call mpp_clock_begin(id_mixed_layer_depth) if(need_data(id_mld, next_time) .or. need_data(id_mld_sq, next_time) & .or. need_data(id_mld_nrho, next_time)) then call 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%neutralrho(isd:ied,jsd:jed,:), & Dens%pressure_at_depth(isd:ied,jsd:jed,:), Time) endif call mpp_clock_end(id_mixed_layer_depth) call mpp_clock_begin(id_mixed_layer_depth_dtheta) if(need_data(id_mld_dtheta, next_time)) then call mixed_layer_depth_dtheta(Thickness, & T_prog(index_temp)%field(isd:ied,jsd:jed,:,tau),& Time) endif call mpp_clock_end(id_mixed_layer_depth_dtheta) call mpp_clock_begin(id_potrho_mixed_layer) call potrho_mixed_layer(Time, Thickness, Dens) call mpp_clock_end(id_potrho_mixed_layer) call mpp_clock_begin(id_compute_subduction) if( need_data(id_subduction, next_time) .or. need_data(id_subduction_dhdt, next_time) & .or. need_data(id_subduction_horz, next_time) .or. need_data(id_subduction_vert, next_time) & .or. need_data(id_subduction_nrho, next_time) .or. need_data(id_subduction_dhdt_nrho, next_time) & .or. need_data(id_subduction_horz_nrho, next_time) .or. need_data(id_subduction_vert_nrho, next_time) & .or. need_data(id_subduction_mld, next_time) .or. compute_subduction_diags) then call compute_subduction(Time, Thickness, Velocity, Adv_vel, Dens, T_prog, & T_prog(index_salt)%field(isd:ied,jsd:jed,:,taup1),& T_prog(index_temp)%field(isd:ied,jsd:jed,:,taup1),& diff_cbt) endif call mpp_clock_end(id_compute_subduction) call mpp_clock_begin(id_compute_tracer_mld) call compute_tracer_mld(Time, Thickness, Dens, T_prog, & T_prog(index_salt)%field(isd:ied,jsd:jed,:,taup1), & T_prog(index_temp)%field(isd:ied,jsd:jed,:,taup1)) call mpp_clock_end(id_compute_tracer_mld) ! compute tracer as function of potential density call mpp_clock_begin(id_diagnose_tracer_on_rho) do n=1,num_prog_tracers if(need_data(id_tracer_on_rho(n), next_time)) then call diagnose_tracer_on_rho(Time, Dens, T_prog(n), n) endif enddo call mpp_clock_end(id_diagnose_tracer_on_rho) ! compute tracer*dz/drho as function of potential density call mpp_clock_begin(id_diagnose_tracer_zrho_on_rho) do n=1,num_prog_tracers if(need_data(id_tracer_zrho_on_rho(n), next_time)) then call diagnose_tracer_zrho_on_rho(Time, Dens, Thickness, T_prog(n), n) endif enddo call mpp_clock_end(id_diagnose_tracer_zrho_on_rho) ! compute depth of isopycnal surfaces and send to diagnostics manager call mpp_clock_begin(id_diagnose_depth_of_potrho) if(need_data(id_depth_of_potrho, next_time)) then call diagnose_depth_of_potrho(Time, Dens, Thickness) endif call mpp_clock_end(id_diagnose_depth_of_potrho) ! compute depth of potential temp surfaces and send to diagnostics manager call mpp_clock_begin(id_diagnose_depth_of_theta) if(need_data(id_depth_of_theta, next_time)) then call diagnose_depth_of_theta(Time, Dens, Thickness, T_prog) endif call mpp_clock_end(id_diagnose_depth_of_theta) call mpp_clock_begin(id_diagnose_mixing) ! quantify levels of effective mixing using sorting algorithm ! method available for tendency==TWO_LEVEL if (tendency==TWO_LEVEL) then if(need_data(id_kappa_sort, next_time) .or. need_data(id_temp_sort,next_time)) then call diagnose_kappa_sort(Time, Thickness, T_prog(index_temp)) endif endif ! quantify levels of mixing using horizontal averaging algorithm if(need_data(id_kappa_simple,next_time)) then call diagnose_kappa_simple(Time, T_prog(index_temp)) endif call mpp_clock_end(id_diagnose_mixing) call mpp_clock_begin(id_tracer_numerical) if (diag_step > 0) then if (mod(Time%itt, diag_step) == 0) then call dzt_min_max(Time, Thickness, 'From ocean_tracer_diag, dzt_min_max information') do n = 1,num_prog_tracers call tracer_min_max(Time, Thickness, T_prog(n)) call tracer_integrals(Time, Thickness, T_prog(n)) enddo call tracer_change(Time, Thickness, T_prog, T_diag, Ext_mode, & pme, melt, runoff, calving) call tracer_land_cell_check (Time, T_prog) endif endif call mpp_clock_end(id_tracer_numerical) call mpp_clock_begin(id_conservation) if (nint(dtts) /= 0) then call mass_conservation (Time, Thickness, Ext_mode, pme, runoff, calving) call tracer_conservation (Time, Thickness, T_prog, T_diag, pme, runoff, calving) endif call mpp_clock_end(id_conservation) ! send total seawater mass to diag_manager call mpp_clock_begin(id_total_mass) if(need_data(id_mass_seawater, next_time)) then call send_total_mass(Time, Thickness) endif call mpp_clock_end(id_total_mass) ! send total seawater volume to diag_manager call mpp_clock_begin(id_total_volume) if(need_data(id_volume_seawater, next_time)) then call send_total_volume(Time, Thickness) endif call mpp_clock_end(id_total_volume) ! global average pressure if(need_data(id_press_ave, next_time)) then call send_global_ave_pressure(Time, Thickness, Dens) endif ! send total tracer to diag_manager call mpp_clock_begin(id_total_tracer) do n=1,num_prog_tracers if(need_data(id_prog_total(n), next_time)) then call send_total_tracer(Time, Thickness, T_prog(n), n) endif if(need_data(id_prog_global_ave(n), next_time)) then call send_global_ave_tracer(Time, Thickness, T_prog(n), n) endif if(need_data(id_prog_surface_ave(n), next_time)) then call send_surface_ave_tracer(Time, Thickness, T_prog(n), n) endif if(need_data(id_prog_surface_area_ave(n), next_time)) then call send_surface_area_ave_tracer(Time, T_prog(n), n) endif enddo call mpp_clock_end(id_total_tracer) end subroutine ocean_tracer_diagnostics ! NAME="ocean_tracer_diagnostics" !####################################################################### ! ! ! ! ! Calculate the mixed layer depth (m), which is defined as the depth ( > 0 ) ! where the buoyancy difference with respect to the surface level is ! equal to buoyancy_crit (m/s^2). ! ! Note that the mixed layer depth is taken with respect to the ocean surface ! at z=eta_t, so the mixed layer depth is always positive. That is, the mld ! is here defined as a thickness of water. ! ! ! subroutine calc_mixed_layer_depth(Thickness, salinity, theta, rho, pressure, & hmxl, smooth_mld_input) type(ocean_thickness_type), intent(in) :: Thickness real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(in) :: rho real, dimension(isd:,jsd:,:), intent(in) :: pressure real, dimension(isd:,jsd:), intent(inout) :: hmxl logical, optional, intent(in) :: smooth_mld_input real, parameter :: epsln=1.0e-20 ! for divisions integer :: i, j, k, km1, kb logical :: smooth_mld_routine if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (calc_mixed_layer_depth): module needs initialization ') endif if (present(smooth_mld_input)) then smooth_mld_routine = smooth_mld_input else smooth_mld_routine = smooth_mld endif hmxl(:,:) = 0.0 wrk1(:,:,:) = 0.0 wrk2(:,:,:) = 0.0 wrk1(:,:,:) = density_delta_sfc( rho(:,:,:), salinity(:,:,:), theta(:,:,:), pressure(:,:,:)) do k=2,nk do j=jsc,jec do i=isc,iec wrk2(i,j,k) = -grav*Grd%tmask(i,j,k)*wrk1(i,j,k-1)/(epsln+rho(i,j,k)) enddo enddo enddo do j=jsc,jec do i=isc,iec kb=Grd%kmt(i,j) if(kb==0) then hmxl(i,j) = 0.0 else hmxl(i,j) = Thickness%depth_zwt(i,j,kb) endif enddo enddo do k=2,nk km1 = k-1 do j=jsc,jec do i=isc,iec kb=Grd%kmt(i,j) if (kb == 0) then hmxl(i,j) = 0.0 else if ( wrk2(i,j,k) >= buoyancy_crit .and. hmxl(i,j)==Thickness%depth_zwt(i,j,kb)) then hmxl(i,j) = Thickness%depth_zt(i,j,km1) & - (Thickness%depth_zt(i,j,km1) - Thickness%depth_zt(i,j,k)) & * (buoyancy_crit-wrk2(i,j,km1)) & / (wrk2(i,j,k) - wrk2(i,j,km1) + epsln) endif hmxl(i,j) = hmxl(i,j) * Grd%tmask(i,j,1) endif enddo enddo enddo ! smooth mld if(smooth_mld_routine) then call mpp_update_domains(hmxl(:,:), Dom%domain2d) hmxl(:,:) = S2D(hmxl(:,:)) endif end subroutine calc_mixed_layer_depth ! NAME="calc_mixed_layer_depth" !####################################################################### ! ! ! ! ! Diagnose mixed layer depth (m). ! Call calc_mixed_layer_depth to determine the mixed layer depth. ! Also compute neutral density at depth of the mixed layer. ! ! subroutine mixed_layer_depth(Thickness, salinity, theta, rho, neutralrho, pressure, Time) type(ocean_thickness_type), intent(in) :: Thickness real, dimension(isd:,jsd:,:), intent(in) :: salinity real, dimension(isd:,jsd:,:), intent(in) :: theta real, dimension(isd:,jsd:,:), intent(in) :: rho real, dimension(isd:,jsd:,:), intent(in) :: neutralrho real, dimension(isd:,jsd:,:), intent(in) :: pressure type(ocean_time_type), intent(in) :: Time integer :: i,j,k,kmt real, dimension(isd:ied,jsd:jed) :: hmxl if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (mixed_layer_depth): module needs initialization ') endif call calc_mixed_layer_depth(Thickness, salinity, theta, rho, pressure, hmxl) call diagnose_2d(Time, Grd, id_mld, hmxl(:,:)) if(id_mld_sq > 0) then call diagnose_2d(Time, Grd, id_mld_sq, hmxl(:,:)**2) endif ! compute neutral density at depth of hmxl if(id_mld_nrho > 0) then wrk1_2d(:,:) = 0.0 do j=jsc,jec do i=isc,iec kmt = max(1,Grd%kmt(i,j)) if(hmxl(i,j) < Thickness%depth_zt(i,j,1)) then wrk1_2d(i,j) = neutralrho(i,j,1) cycle endif if(hmxl(i,j) > Thickness%depth_zt(i,j,kmt)) then wrk1_2d(i,j) = neutralrho(i,j,kmt) cycle endif do k=2,nk if(hmxl(i,j) > Thickness%depth_zt(i,j,k-1) .and. & hmxl(i,j) <= Thickness%depth_zt(i,j,k)) then wrk1_2d(i,j) = ( neutralrho(i,j,k-1)*(Thickness%depth_zt(i,j,k)-hmxl(i,j)) & +neutralrho(i,j,k) *(hmxl(i,j)-Thickness%depth_zt(i,j,k-1))& )/Thickness%dzwt(i,j,k-1) endif enddo enddo enddo call diagnose_2d(Time, Grd, id_mld_nrho, wrk1_2d(:,:)) endif end subroutine mixed_layer_depth ! NAME="mixed_layer_depth" !####################################################################### ! ! ! ! ! Calculate the depth required to reach a temperature that is ! dtheta cooler than the surface temperature. ! ! Note: ! 1/ mixed_layer_depth_dtheta is taken with respect to the ocean surface ! at z=eta_t. ! ! 2/ mixed_layer_depth_dtheta is no greater than the ocean depth + eta_t. ! ! Coded March 2010 by Stephen.Griffies ! ! ! subroutine mixed_layer_depth_dtheta(Thickness, theta, Time) type(ocean_thickness_type), intent(in) :: Thickness real, dimension(isd:,jsd:,:), intent(in) :: theta type(ocean_time_type), intent(in) :: Time real :: epsln=1.0e-20 ! for divisions real :: dtheta integer :: i, j, k, km1, kb if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (mixed_layer_depth_dtheta): module needs initialization ') endif wrk1(:,:,:) = 0.0 wrk1_2d(:,:) = 0.0 do j=jsc,jec do i=isc,iec kb=Grd%kmt(i,j) if(kb==0) then wrk1_2d(i,j) = 0.0 else wrk1_2d(i,j) = Thickness%depth_zwt(i,j,kb) endif enddo enddo do k=2,nk km1=k-1 do j=jsc,jec do i=isc,iec kb=Grd%kmt(i,j) if (kb == 0) then wrk1_2d(i,j) = 0.0 else dtheta = theta(i,j,1) - theta(i,j,k) if (dtheta >= dtheta_crit .and. wrk1_2d(i,j)==Thickness%depth_zwt(i,j,kb)) then wrk1_2d(i,j) = Thickness%depth_zt(i,j,km1) & + (dtheta_crit/(dtheta+epsln)) & *(Thickness%depth_zt(i,j,k) - Thickness%depth_zt(i,j,km1)) endif wrk1_2d(i,j) = wrk1_2d(i,j) * Grd%tmask(i,j,1) endif enddo enddo enddo call diagnose_2d(Time, Grd, id_mld_dtheta, wrk1_2d(:,:)) end subroutine mixed_layer_depth_dtheta ! NAME="mixed_layer_depth_dtheta" !####################################################################### ! ! ! ! Diagnose subduction rate (kg/sec) based on kinematic method to compute ! mass transport through base of mixed layer. ! ! Some approximations made for convenience: ! 1/ use velocity at time tau ! use tracer at time tau ! use thickness at time taup1, since all pieces of Thickness ! have already been updated to taup1. ! ! 2/ horizontally interpolate B-grid to C-grid velocity components, ! but then vertically interpolate using mld computed at T-points. ! ! Stephen.Griffies ! March 2012 ! Updated for tracer transport July 2013 ! ! subroutine compute_subduction(Time, Thickness, Velocity, Adv_vel, Dens, T_prog, salinity, theta, diff_cbt) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness ! taup1 values type(ocean_velocity_type), intent(in) :: Velocity type(ocean_adv_vel_type), intent(in) :: Adv_vel ! tau values type(ocean_density_type), intent(inout) :: Dens type(ocean_prog_tracer_type), intent(in) :: T_prog(:) ! for other than temp/saln tracers real, dimension(isd:,jsd:,:), intent(in) :: salinity ! taup1 value real, dimension(isd:,jsd:,:), intent(in) :: theta ! taup1 value real, dimension(isd:,jsd:,:,:), intent(in) :: diff_cbt integer :: i,j,k,kmt,n,nmix integer :: tau, taup1 real :: denominator, denominator_r real :: W1, W2 real :: dhdx, dhdy, dhdt real, dimension(isd:ied,jsd:jed) :: mld_tau real, dimension(isd:ied,jsd:jed) :: mld_taup1 real, dimension(isd:ied,jsd:jed) :: rho_mld real, dimension(isd:ied,jsd:jed) :: subduction real, dimension(isd:ied,jsd:jed) :: subduction_dhdt real, dimension(isd:ied,jsd:jed) :: subduction_horz real, dimension(isd:ied,jsd:jed) :: subduction_vert real, dimension(isd:ied,jsd:jed,3) :: velocity_mld if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (compute_subduction): module needs initialization') endif tau = Time%tau taup1 = Time%taup1 ! compute mld_taup1 ! note Thickness is filled with taup1 thickness fields. call calc_mixed_layer_depth(Thickness, salinity, theta, & Dens%rho(isd:ied,jsd:jed,:,taup1), & Dens%pressure_at_depth(isd:ied,jsd:jed,:), & mld_taup1(:,:), & smooth_mld_for_subduction) ! update Dens%mld_subduction do j=jsc,jec do i=isc,iec mld_tau(i,j) = Dens%mld_subduction(i,j) Dens%mld_subduction(i,j) = mld_taup1(i,j) enddo enddo ! interpolate B-grid horizontal velocity to C-grid wrk1_v(:,:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec denominator = Grd%dus(i,j) + Grd%dun(i,j-1) + epsln wrk1_v(i,j,k,1) = Grd%tmask(i+1,j,k)*(Velocity%u(i,j,k,1,tau) *Grd%dun(i,j-1) & + Velocity%u(i,j-1,k,1,tau)*Grd%dus(i,j))/denominator denominator = Grd%duw(i,j) + Grd%due(i-1,j) + epsln wrk1_v(i,j,k,2) = Grd%tmask(i,j+1,k)*(Velocity%u(i-1,j,k,2,tau)*Grd%duw(i,j) & + Velocity%u(i,j,k,2,tau) *Grd%due(i-1,j))/denominator enddo enddo enddo ! vertical velocity component at T-cell wrk1(:,:,:) = 0.0 if(vert_coordinate_class==DEPTH_BASED) then do k=1,nk do j=jsc,jec do i=isc,iec wrk1(i,j,k) = rho0r*Adv_vel%wrho_bt(i,j,k) enddo enddo enddo else do k=1,nk do j=jsc,jec do i=isc,iec if(Grd%tmask(i,j,k) > 0.0) then wrk1(i,j,k) = Adv_vel%wrho_bt(i,j,k)/Dens%rho(i,j,k,tau) endif enddo enddo enddo endif ! vertically interpolate velocity components to mld_tau; ! ignore horizontal offset between C-grid velocity and mld. velocity_mld(:,:,:) = 0.0 rho_mld(:,:) = 0.0 ! shallow mld_tau do j=jsc,jec do i=isc,iec if(mld_tau(i,j) <= Thickness%depth_zt(i,j,1)) then velocity_mld(i,j,1) = wrk1_v(i,j,1,1) velocity_mld(i,j,2) = wrk1_v(i,j,1,2) velocity_mld(i,j,3) = wrk1(i,j,1) rho_mld(i,j) = Dens%rho(i,j,1,tau) endif enddo enddo ! intermediate mld_tau do j=jsc,jec do i=isc,iec kloop: do k=1,nk-1 if(Thickness%depth_zt(i,j,k) < mld_tau(i,j) .and. mld_tau(i,j) <= Thickness%depth_zt(i,j,k+1)) then if(Grd%tmask(i,j,k+1) > 0) then W1= mld_tau(i,j) - Thickness%depth_zt(i,j,k) W2= Thickness%depth_zt(i,j,k+1) - mld_tau(i,j) denominator_r = 1.0/(W1 + W2 + epsln) velocity_mld(i,j,1) = (wrk1_v(i,j,k,1)*W2 + wrk1_v(i,j,k+1,1)*W1)*denominator_r velocity_mld(i,j,2) = (wrk1_v(i,j,k,2)*W2 + wrk1_v(i,j,k+1,2)*W1)*denominator_r velocity_mld(i,j,3) = (wrk1(i,j,k) *W2 + wrk1(i,j,k+1) *W1)*denominator_r rho_mld(i,j) = (Dens%rho(i,j,k,tau)*W2 + Dens%rho(i,j,k+1,tau)*W1)*denominator_r exit kloop endif endif enddo kloop enddo enddo ! deep mld_tau do j=jsc,jec do i=isc,iec kmt = Grd%kmt(i,j) if(kmt > 0) then if(mld_tau(i,j) > Thickness%depth_zt(i,j,kmt)) then velocity_mld(i,j,1) = wrk1_v(i,j,kmt,1) velocity_mld(i,j,2) = wrk1_v(i,j,kmt,2) velocity_mld(i,j,3) = wrk1(i,j,kmt) rho_mld(i,j) = Dens%rho(i,j,kmt,tau) endif endif enddo enddo ! over-write rho_mld for Boussinesq case if(vert_coordinate_class==DEPTH_BASED) then do j=jsc,jec do i=isc,iec rho_mld(i,j) = rho0*Grd%tmask(i,j,1) enddo enddo endif ! save rho*dat wrk1_2d(:,:) = 0.0 do j=jsc,jec do i=isc,iec wrk1_2d(i,j) = rho_mld(i,j)*Grd%dat(i,j) enddo enddo ! compute contributions to subduction call mpp_update_domains(mld_tau(:,:), Dom%domain2d) subduction_dhdt(:,:) = 0.0 subduction_horz(:,:) = 0.0 subduction_vert(:,:) = 0.0 subduction(:,:) = 0.0 do j=jsc,jec do i=isc,iec dhdx = (mld_tau(i+1,j)-mld_tau(i,j))*Grd%dxter(i,j) dhdy = (mld_tau(i,j+1)-mld_tau(i,j))*Grd%dytnr(i,j) dhdt = (mld_taup1(i,j) - mld_tau(i,j))*dtimer subduction_dhdt(i,j) = -wrk1_2d(i,j)*dhdt subduction_horz(i,j) = -wrk1_2d(i,j)*(dhdx*velocity_mld(i,j,1)+dhdy*velocity_mld(i,j,2)) subduction_vert(i,j) = -wrk1_2d(i,j)*velocity_mld(i,j,3) subduction(i,j) = subduction_dhdt(i,j) + subduction_horz(i,j) + subduction_vert(i,j) enddo enddo call diagnose_2d(Time, Grd, id_subduction_dhdt, subduction_dhdt(:,:)) call diagnose_2d(Time, Grd, id_subduction_horz, subduction_horz(:,:)) call diagnose_2d(Time, Grd, id_subduction_vert, subduction_vert(:,:)) call diagnose_2d(Time, Grd, id_subduction, subduction(:,:)) call diagnose_2d(Time, Grd, id_subduction_mld, mld_tau(:,:)) if(id_subduction_dhdt_nrho > 0) then wrk1(:,:,:) = 0.0 k=1 do j=jsc,jec do i=isc,iec wrk1(i,j,k) = subduction_dhdt(i,j) enddo enddo call diagnose_3d_rho(Time, Dens, id_subduction_dhdt_nrho, wrk1) endif if(id_subduction_horz_nrho > 0) then wrk1(:,:,:) = 0.0 k=1 do j=jsc,jec do i=isc,iec wrk1(i,j,k) = subduction_horz(i,j) enddo enddo call diagnose_3d_rho(Time, Dens, id_subduction_horz_nrho, wrk1) endif if(id_subduction_vert_nrho > 0) then wrk1(:,:,:) = 0.0 k=1 do j=jsc,jec do i=isc,iec wrk1(i,j,k) = subduction_vert(i,j) enddo enddo call diagnose_3d_rho(Time, Dens, id_subduction_vert_nrho, wrk1) endif if(id_subduction_nrho > 0) then wrk1(:,:,:) = 0.0 k=1 do j=jsc,jec do i=isc,iec wrk1(i,j,k) = subduction(i,j) enddo enddo call diagnose_3d_rho(Time, Dens, id_subduction_nrho, wrk1) endif !!!!!!! tracer diagnostic calculations !!!!!!! do n=1,num_prog_tracers if (n==index_salt) then nmix=2 else nmix=1 endif ! compute downgradient vertical diffusive flux (assume time-implicit calculation) wrk1(:,:,:) = 0.0 if(id_subduction_mld_zflux_diff(n) > 0) then do k=1,nk-1 do j=jsc,jec do i=isc,iec wrk1(i,j,k) = -1.0*rho0*Grd%tmask(i,j,k+1)*diff_cbt(i,j,k,nmix) & *(T_prog(n)%field(i,j,k,taup1)-T_prog(n)%field(i,j,k+1,taup1))& /Thickness%dzwt(i,j,k) enddo enddo enddo endif if( id_subduction_tracer(n) >0 .or. id_subduction_dhdt_tracer(n) >0 .or. id_subduction_horz_tracer(n)>0 & .or. id_subduction_vert_tracer(n)>0 .or. id_subduction_mld_zflux_diff(n)>0) then ! interpolate tracer concentration and vertical diffusive flux to base of mld wrk1_2d(:,:) = 0.0 ! interpolated tracer wrk2_2d(:,:) = 0.0 ! interpolated zflux ! shallow mld_tau do j=jsc,jec do i=isc,iec if(mld_tau(i,j) <= Thickness%depth_zt(i,j,1)) then wrk1_2d(i,j) = T_prog(n)%field(i,j,k,tau) wrk2_2d(i,j) = wrk1(i,j,k) endif enddo enddo ! intermediate mld_tau do j=jsc,jec do i=isc,iec kloopB: do k=1,nk-1 if(Thickness%depth_zt(i,j,k) < mld_tau(i,j) .and. mld_tau(i,j) <= Thickness%depth_zt(i,j,k+1)) then if(Grd%tmask(i,j,k+1) > 0) then W1= mld_tau(i,j) - Thickness%depth_zt(i,j,k) W2= Thickness%depth_zt(i,j,k+1) - mld_tau(i,j) denominator_r = 1.0/(W1 + W2 + epsln) wrk1_2d(i,j) = (T_prog(n)%field(i,j,k,tau)*W2 + T_prog(n)%field(i,j,k+1,tau)*W1)*denominator_r wrk2_2d(i,j) = (wrk1(i,j,k)*W2 + wrk1(i,j,k+1)*W1)*denominator_r exit kloopB endif endif enddo kloopB enddo enddo ! deep mld_tau do j=jsc,jec do i=isc,iec kmt = Grd%kmt(i,j) if(kmt > 0) then if(mld_tau(i,j) > Thickness%depth_zt(i,j,kmt)) then wrk1_2d(i,j) = T_prog(n)%field(i,j,kmt,tau) wrk2_2d(i,j) = wrk1(i,j,kmt) endif endif enddo enddo endif ! write diagnostics wrk3_2d(:,:) = 0.0 if(id_subduction_tracer(n) > 0) then wrk3_2d(:,:) = subduction(:,:)*wrk1_2d(:,:) used = send_data (id_subduction_tracer(n), wrk3_2d(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif if(id_subduction_dhdt_tracer(n) > 0) then wrk3_2d(:,:) = subduction_dhdt(:,:)*wrk1_2d(:,:) used = send_data (id_subduction_dhdt_tracer(n), wrk3_2d(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif if(id_subduction_horz_tracer(n) > 0) then wrk3_2d(:,:) = subduction_horz(:,:)*wrk1_2d(:,:) used = send_data (id_subduction_horz_tracer(n), wrk3_2d(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif if(id_subduction_vert_tracer(n) > 0) then wrk3_2d(:,:) = subduction_vert(:,:)*wrk1_2d(:,:) used = send_data (id_subduction_vert_tracer(n), wrk3_2d(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif if(id_subduction_mld_zflux_diff(n) > 0) then wrk3_2d(:,:) = wrk2_2d(:,:)*Grd%dat(:,:) used = send_data (id_subduction_mld_zflux_diff(n), wrk3_2d(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif enddo ! end of n-loop end subroutine compute_subduction ! NAME="compute_subduction" !####################################################################### ! ! ! ! ! Diagnose the vertically integrated tracer within the mixed layer. ! ! Work with taup1 values, since Thickness type is filled with taup1 ! thickness fields. ! ! Stephen.Griffies ! July 2013 ! ! subroutine compute_tracer_mld(Time, Thickness, Dens, T_prog, salinity, theta) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness ! taup1 values type(ocean_density_type), intent(inout) :: Dens type(ocean_prog_tracer_type), intent(in) :: T_prog(:) real, dimension(isd:,jsd:,:), intent(in) :: salinity ! taup1 value real, dimension(isd:,jsd:,:), intent(in) :: theta ! taup1 value 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, & '==>Error from ocean_tracer_diag_mod (compute_tracer_mld): module needs initialization') endif if(.not. compute_tracer_mld_diags) return taup1 = Time%taup1 ! compute mld using taup1 values. ! Note the Thickness type is filled with taup1 thickness fields. call calc_mixed_layer_depth(Thickness, salinity, theta, & Dens%rho(isd:ied,jsd:jed,:,taup1), & Dens%pressure_at_depth(isd:ied,jsd:jed,:), & mld(:,:), smooth_mld_for_subduction) ! determine column masking function 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 ! k>1 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 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 ! compute the vertically integrated tracer within the mixed layer do n=1,num_prog_tracers if(id_tracer_mld(n) > 0) then wrk1_2d(:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec wrk1_2d(i,j) = wrk1_2d(i,j) & + wrk1(i,j,k)*Thickness%rho_dzt(i,j,k,taup1)*T_prog(n)%field(i,j,k,taup1) enddo enddo enddo used = send_data (id_tracer_mld(n), wrk1_2d(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif enddo end subroutine compute_tracer_mld ! NAME="compute_tracer_mld" !####################################################################### ! ! ! ! ! Compute change in tracer over a time step and difference between ! this change and the boundary forcing. ! ! This routine is very useful for detecting bugs in tracer routines. ! ! subroutine tracer_change (Time, Thickness, T_prog, T_diag, Ext_mode, & pme, melt, runoff, calving) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: T_prog(:) type(ocean_diag_tracer_type), intent(in) :: T_diag(:) type(ocean_external_mode_type), intent(in) :: Ext_mode real, dimension(isd:,jsd:), intent(in) :: pme real, dimension(isd:,jsd:), intent(in) :: melt real, dimension(isd:,jsd:), intent(in) :: runoff real, dimension(isd:,jsd:), intent(in) :: calving real, dimension(isd:ied,jsd:jed) :: tmp real, dimension(isd:ied,jsd:jed) :: tmpE real, dimension(isd:ied,jsd:jed) :: tmpL real, dimension(isd:ied,jsd:jed) :: area_k real :: mass_error, tracer_error, flux_error real :: tracer_input_stf, tracer_input_btf, tracer_input_pme real :: tracer_input_otf real :: tracer_input_eta_smooth, tracer_input_pbot_smooth real :: tracer_input_runoff, tracer_input_calving, tracer_input_total real :: pme_input, runoff_input, calving_input, melt_input, obc_input real :: eta_t_max, eta_t_min, eta_t_avg real :: anompbot_max, anompbot_min real :: frazil_heat real :: tracer_th_tend, tracer_src real :: ttracer real :: tracer_L, tracer_E real :: dmasstracer real :: dmasstracer_E, dmasstracer_L, dmasstracer_L2 real :: mass_total, mass_source, mass_change real :: mass_E, mass_L, mass_T real :: mass_changeE, mass_changeL, mass_changeL2 real :: mass_taup1, mass_taup1_r, mass_taup1_keq1 real :: mass_eta_smooth, mass_pbot_smooth real :: conversion, conv_blob integer :: i,j,k,n integer :: taum1, tau, taup1 integer :: stdoutunit stdoutunit=stdout() call mpp_clock_begin(id_tracer_change) if (.not.module_is_initialized) then call mpp_error(FATAL,& '==>Error from ocean_tracer_diag_mod (tracer_change): module needs initialization ') endif if(size(T_prog,1) /= num_prog_tracers) then call mpp_error(FATAL,& '==>Error from ocean_tracer_diag_mod (tracer_change): passing T_prog with wrong dimensions') endif write(stdoutunit,*) ' ' write (stdoutunit,'(//50x,a)') ' Mass and tracer change summary (over one timestep):' call write_timestamp(Time%model_time) write(stdoutunit,*) ' ' ! compute total mass and changes over a time step pme_input = 0.0 ! total mass of evap-precip input to ocean over time step (kg) melt_input = 0.0 ! total mass of ice melt input to ocean over time step (kg) runoff_input = 0.0 ! total mass of runoff water input to ocean over time step (kg) calving_input = 0.0 ! total mass of calving water input to ocean over time step (kg) obc_input = 0.0 ! total mass of water input through open boundaries to ocean over time step (kg) eta_t_avg = 0.0 ! area average of eta_t (m) eta_t_max = 0.0 ! global max of eta_t (m) eta_t_min = 0.0 ! global min of eta_t (m) anompbot_max = 0.0 ! global max of pbot_t-pbot0 (dbar) anompbot_min = 0.0 ! global min of pbot_t-pbot0 (dbar) mass_total = 0.0 ! mass of ocean on tracer cells (kg) mass_taup1_keq1 = 0.0 ! mass of ocean in the k=1 cell (kg) mass_eta_smooth = 0.0 ! mass from eta_t surface filtering (kg) mass_pbot_smooth = 0.0 ! mass from pbot_t filtering (kg) mass_change = 0.0 ! change in mass over time step (kg) mass_source = 0.0 ! mass from sources (kg) mass_L = 0.0 ! mass of L system mass_E = 0.0 ! mass of E system mass_T = 0.0 mass_changeL = 0.0 ! change in mass of L system mass_changeE = 0.0 ! change in mass of E system mass_changeL2 = 0.0 ! mass transfer from L to E system ! note: if all is well, mass_changeL = mass_changeL2 conv_blob = 0.0 ! tendency due to convergence of blobs (should be close to zero) taum1 = Time%taum1 tau = Time%tau taup1 = Time%taup1 ! inverse mass of ocean water at time taup1 mass_taup1 = total_mass(Thickness, taup1) mass_taup1_r = 1/(mass_taup1+epsln) ! mass of ocean water at time taup1 and k-level=1 mass_taup1_keq1 = klevel_total_mass(Thickness, taup1, 1) ! mass of ocean water at time tau tmpE(:,:) = 0.0 do k=1,nk tmpE(:,:) = tmpE(:,:) + Grd%tmask(:,:,k)*area_t(:,:)*Thickness%rho_dzt(:,:,k,tau) enddo mass_E = mpp_global_sum(Dom%domain2d, tmpE(:,:), global_sum_flag) mass_total = mass_E ! mass of ocean water at time tau, no with blobs if(use_blobs) then tmpL(:,:) = 0.0 tmp(:,:) = 0.0 do k=1,nk tmpL(:,:) = tmpL(:,:) + Grd%tmask(:,:,k)*area_t(:,:)*Thickness%rho_dztL(:,:,k,tau) tmp(:,:) = tmp(:,:) + Grd%tmask(:,:,k)*area_t(:,:)*Thickness%rho_dztT(:,:,k,tau) enddo mass_L = mpp_global_sum(Dom%domain2d, tmpL(:,:), global_sum_flag) mass_T = mpp_global_sum(Dom%domain2d, tmp(:,:), global_sum_flag) mass_total = mass_L + mass_E endif ! mass of ice melt input to the ocean tmp(:,:) = dtime*melt(:,:)*area_t(:,:) melt_input = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! mass of pme input to the ocean. subtract melt_input ! since ice melt is already contained in pme. tmp(:,:) = dtime*pme(:,:)*area_t(:,:) pme_input = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) - melt_input ! mass of river runoff water input to the ocean tmp(:,:) = dtime*runoff(:,:)*area_t(:,:) runoff_input = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! mass of calving land ice input to the ocean tmp(:,:) = dtime*calving(:,:)*area_t(:,:) calving_input = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! area average surface height at time tau tmp(:,:) = area_t(:,:)*Ext_mode%eta_t(:,:,tau) eta_t_avg = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag)/Grd%tcellsurf ! max and min of eta_t eta_t_max = mpp_global_max(Dom%domain2d,Ext_mode%eta_t(:,:,tau)) eta_t_min = mpp_global_min(Dom%domain2d,Ext_mode%eta_t(:,:,tau)) ! max and min of pbot_t-pbot0 anompbot_max = c2dbars & *mpp_global_max(Dom%domain2d,Ext_mode%pbot_t(:,:,tau)-Thickness%pbot0(:,:)) anompbot_min = c2dbars & *mpp_global_min(Dom%domain2d,Ext_mode%pbot_t(:,:,tau)-Thickness%pbot0(:,:)) ! change in mass over a time step tmpE(:,:) = 0.0 do k=1,nk tmpE(:,:) = tmpE(:,:) + Grd%tmask(:,:,k)*area_t(:,:) & *(Thickness%rho_dzt(:,:,k,taup1)-Thickness%rho_dzt(:,:,k,taum1)) enddo mass_changeE = mpp_global_sum(Dom%domain2d,tmpE(:,:), global_sum_flag) mass_change = mass_changeE ! mass change over a time step with blobs if (use_blobs) then tmpL(:,:) = 0.0 do k=1,nk tmpL(:,:) = tmpL(:,:) + Grd%tmask(:,:,k)*area_t(:,:) & *(Thickness%rho_dztL(:,:,k,taup1)-Thickness%rho_dztL(:,:,k,taum1)) enddo mass_changeL = mpp_global_sum(Dom%domain2d,tmpL(:,:), global_sum_flag) mass_change = mass_changeE + mass_changeL tmp(:,:) = dtime*Grd%tmask(:,:,1)*Grd%dat(:,:)*Thickness%blob_source(:,:) mass_changeL2 = mass_changeL2 + mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = dtime*area_t(:,:)*Ext_mode%conv_blob(:,:) conv_blob = mpp_global_sum(Dom%domain2d, tmp(:,:), global_sum_flag) endif ! mass input to ocean from sources tmp(:,:) = 0.0 do k=1,nk tmp(:,:) = tmp(:,:) + dtime*area_t(:,:)*Thickness%mass_source(:,:,k) enddo mass_source = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! mass input due to smoothing of the surface height (should be roundoff). ! note that this contribution is already in mass_source tmp(:,:) = dtime*area_t(:,:)*Ext_mode%eta_smooth(:,:) mass_eta_smooth = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! mass input due to smoothing of the bottom pressure (should be roundoff). ! note that this contribution is already in mass_source tmp(:,:) = dtime*area_t(:,:)*Ext_mode%pbot_smooth(:,:) mass_pbot_smooth = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! compute mass passing across open boundaries if (have_obc) then call ocean_obc_mass_flux(Time, Ext_mode, tmp) tmp(:,:) = tmp(:,:)*dtime obc_input = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) endif mass_error = mass_change-pme_input-melt_input-runoff_input-calving_input-obc_input-mass_source write (stdoutunit,'(/a)') ' ----Single time step diagnostics for volume and mass----' write (stdoutunit,'(a,es24.17,a)') ' Global max (pbot_t-pbot0) on tracer cells (tau) = ',& anompbot_max,' dbar' write (stdoutunit,'(a,es24.17,a)') ' Global min (pbot_t-pbot0) on tracer cells (tau) = ',& anompbot_min,' dbar' write (stdoutunit,'(a,es24.17,a)') ' Global max surface height on tracer cells (tau) = ',& eta_t_max,' m' write (stdoutunit,'(a,es24.17,a)') ' Global min surface height on tracer cells (tau) = ',& eta_t_min,' m' write (stdoutunit,'(a,es24.17,a)') ' Area average surface height on tracer cells (tau) = ',& eta_t_avg,' m' write (stdoutunit,'(a,es24.17,a)') ' Surface area of k=1 tracer cells = ',& Grd%tcellsurf,' m^2' write (stdoutunit,'(a,es24.17,a)') ' Area integral of surface height on tracer cells = ',& eta_t_avg*Grd%tcellsurf,' m^3' if (use_blobs) then write (stdoutunit,'(a,es24.17,a)') ' L system mass change (taup1 - taum1) = ',& mass_changeL,' kg' write (stdoutunit,'(a,es24.17,a)') ' E system mass change (taup1 - taum1) = ',& mass_changeE,' kg' write (stdoutunit,'(a,es24.17,a)') ' Transfer of mass from L to E (taup1) = ',& mass_changeL2,' kg' endif write (stdoutunit,'(a,es24.17,a)') ' Total Mass change = ',& mass_change,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of ocean tracer cells (tau) = ',& mass_total,' kg' if (use_blobs) then write (stdoutunit,'(a,es24.17,a)') ' Mass of L system ocean tracer cells (tau) = ',& mass_L,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of E system ocean tracer cells (tau) = ',& mass_E,' kg' write (stdoutunit,'(a,es24.17,a)') ' Total Mass of ocean tracer cells (tau) = ',& mass_total,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mismatch between T Mass and L + E mass (T-[L+E]) = ',& (mass_T - (mass_total)),' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass from blob convergence (should be small) = ',& conv_blob,' kg' endif write (stdoutunit,'(a,es24.17,a)') ' Mass from sources = ',& mass_source,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass from eta smoothing = ',& mass_eta_smooth,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass from pbot smoothing = ',& mass_pbot_smooth,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of precip-evap input = ',& pme_input,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of river runoff liquid water input = ',& runoff_input,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of sea ice melt input = ',& melt_input,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of calving land ice input = ',& calving_input,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass flux through open boundaries = ',& obc_input,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mismatch between mass change and water input = ',& mass_error,' kg' ! compute total tracer and changes over a time step do n=1,num_prog_tracers dmasstracer = 0.0 ! (rho*vol)(taup1)*tracer(taup1)-(rho*vol)(taum1)*tracer(taum1)) tracer_input_stf = 0.0 ! tracer quantity input to ocean via stf (>0 means tracer entering ocean) tracer_input_btf = 0.0 ! tracer quantity input to ocean via btf (>0 means tracer leaving ocean) tracer_input_otf = 0.0 ! tracer quantity input through open boundaries (>0 means tracer entering ocean) tracer_input_pme = 0.0 ! tracer quantity input to ocean via pme flux tracer_input_runoff = 0.0 ! tracer quantity input to ocean via river runoff tracer_input_calving= 0.0 ! tracer quantity input to ocean via calving land ice tracer_input_total = 0.0 ! total tracer quantity input to ocean frazil_heat = 0.0 ! heat (Joule) from frazil formation tracer_th_tend = 0.0 ! contribution from tracer tendencies including ! neutral diffusion, sigma diffusion, overflow, xlandmix, ! sponges, rivers, and other more traditional sources. tracer_src = 0.0 ! sources that have dimensions tracer concentration per time tracer_input_eta_smooth = 0.0 ! tracer input at surface due to smoothing eta_t (should be roundoff) tracer_input_pbot_smooth= 0.0 ! tracer input bottom due to smoothing pbot_t (should be roundoff) tracer_L = 0.0 ! tracer quantity in the L system tracer_E = 0.0 ! tracer quantity in the E system dmasstracer_L = 0.0 ! change in tracer quantity of the L system dmasstracer_L2 = 0.0 ! tracer quantity transferred from the L to E system dmasstracer_E = 0.0 ! change in tracer quantity of the E system conversion = T_prog(n)%conversion tmp(:,:) = area_t(:,:)*conversion*dtime*T_prog(n)%stf(:,:) tracer_input_stf = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = area_t(:,:)*conversion*dtime*T_prog(n)%btf(:,:) tracer_input_btf = -mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = area_t(:,:)*conversion*dtime*pme(:,:)*T_prog(n)%tpme(:,:) tracer_input_pme = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! runoff can be negative, as for the "geoengineering" trick of siphoning ! Arabian Sea water to the Caspian Sea, to keep the Caspian from evaporating ! away during a coupled climate simulation. As negative runoff is ignored ! in the rivermix module (rightly so, as we do not wish to insert tracer ! from negative runoff into the ocean), we also wish to ignore negative ! runoff for the diagnostics here, in order to reflect what the prognostic ! equations are doing. tmp(:,:) = 0.0 do j=jsc,jec do i=isc,iec if(runoff(i,j) > 0) then tmp(i,j) = area_t(i,j)*conversion*dtime*runoff(i,j)*T_prog(n)%trunoff(i,j) endif enddo enddo tracer_input_runoff = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = area_t(:,:)*conversion*dtime*calving(:,:)*T_prog(n)%tcalving(:,:) tracer_input_calving = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = area_t(:,:)*conversion*dtime*T_prog(n)%eta_smooth(:,:) tracer_input_eta_smooth = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = area_t(:,:)*conversion*dtime*T_prog(n)%pbot_smooth(:,:) tracer_input_pbot_smooth = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) ! ocean_obc_tracer_flux returns vertically integrated horizontal tracer flux ! at last internal points and zero otherwise if(have_obc) then call ocean_obc_tracer_flux(T_prog(n), tmp) tmp(:,:) = tmp(:,:)*conversion*dtime tracer_input_otf = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) endif ! river(=calving+runoff) is added via T_prog%th_tendency inside ocean_rivermix_mod tracer_input_total = tracer_input_stf + tracer_input_btf + tracer_input_otf & + tracer_input_runoff + tracer_input_calving + tracer_input_pme & + tracer_input_eta_smooth + tracer_input_pbot_smooth frazil_heat = 0.0 if(T_prog(n)%name =='temp' .and. index_frazil > 0) then tmp(:,:) = 0.0 do k=1,nk tmp(:,:) = tmp(:,:) + T_diag(index_frazil)%field(:,:,k)*area_t(:,:)*Grd%tmask(:,:,k)/frazil_factor enddo frazil_heat = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tracer_input_total = tracer_input_total + frazil_heat endif ! tracer changes and total source tmp(:,:) = 0.0 tmpE(:,:) = 0.0 wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 do k=1,nk area_k(:,:) = Grd%dat(:,:)*Grd%tmask(:,:,k) tmpE(:,:) = tmpE(:,:) + area_k(:,:)*conversion* & (Thickness%rho_dzt(:,:,k,taup1)*T_prog(n)%field(:,:,k,taup1) & -Thickness%rho_dzt(:,:,k,taum1)*T_prog(n)%field(:,:,k,taum1)) wrk1_2d(:,:) = wrk1_2d(:,:) + area_k(:,:)*conversion*dtime*T_prog(n)%th_tendency(:,:,k) wrk2_2d(:,:) = wrk2_2d(:,:) + area_k(:,:)*conversion*dtime*T_prog(n)%source(:,:,k)*Thickness%rho_dzt(:,:,k,taup1) enddo if(have_obc) then tmpE(:,:) = tmpE(:,:)*Grd%obc_tmask(:,:) wrk1_2d(:,:) = wrk1_2d(:,:)*Grd%obc_tmask(:,:) wrk2_2d(:,:) = wrk2_2d(:,:)*Grd%obc_tmask(:,:) endif tracer_th_tend = mpp_global_sum(Dom%domain2d, wrk1_2d(:,:), global_sum_flag) tracer_src = mpp_global_sum(Dom%domain2d, wrk2_2d(:,:), global_sum_flag) dmasstracer_E = mpp_global_sum(Dom%domain2d, tmpE(:,:) , global_sum_flag) dmasstracer = dmasstracer_E ! tracer changes and total source, now with blobs if(use_blobs) then tmpL(:,:) = 0.0 tmp(:,:) = 0.0 do k=1,nk area_k(:,:) = Grd%dat(:,:)*Grd%tmask(:,:,k) tmpL(:,:) = tmpL(:,:) + conversion*(T_prog(n)%sum_blob(:,:,k,taup1) - T_prog(n)%sum_blob(:,:,k,taum1)) tmp(:,:) = tmp(:,:) + dtime*conversion*area_k(:,:)*(T_prog(n)%tend_blob(:,:,k)) enddo if(have_obc) then tmpL(:,:) = tmpL(:,:)*Grd%obc_tmask(:,:) tmp(:,:) = tmp(:,:)*Grd%obc_tmask(:,:) endif dmasstracer_L = dmasstracer_L + mpp_global_sum(Dom%domain2d, tmpL(:,:), global_sum_flag) dmasstracer_L2 = dmasstracer_L2 + mpp_global_sum(Dom%domain2d, tmp(:,:), global_sum_flag) dmasstracer = dmasstracer_L + dmasstracer_E endif ! for cleaner diagnostics, remove contributions from runoff, calving, pme, and smoother, ! each of which are added to T_prog%th_tendency inside of ocean_tracer.F90. tracer_th_tend = tracer_th_tend -tracer_input_runoff -tracer_input_calving -tracer_input_pme & -tracer_input_eta_smooth -tracer_input_pbot_smooth ! for aidif=0.0, stf is included in T_prog%th_tendency through vertical diffusion. ! for aidif=1.0, it is not in th_tendency; rather, it is included through invtri. if(aidif==0.0) then tracer_th_tend = tracer_th_tend - tracer_input_stf endif ! tracer in E system tmpE(:,:) = 0.0 do k=1,nk tmpE(:,:) = tmpE(:,:) + conversion*area_t(:,:)*Thickness%rho_dzt(:,:,k,tau)*T_prog(n)%field(:,:,k,tau) enddo tracer_E = mpp_global_sum(Dom%domain2d, tmpE(:,:), global_sum_flag) ttracer = tracer_E ! tracer in L system if (use_blobs) then tmpL(:,:) = 0.0 do k=1,nk tmpL(:,:) = tmpL(:,:) + conversion*Grd%tmask(:,:,k)*(T_prog(n)%sum_blob(:,:,k,tau)) enddo tracer_L = mpp_global_sum(Dom%domain2d, tmpL(:,:), global_sum_flag) ttracer = tracer_L + tracer_E endif tracer_error = dmasstracer-tracer_input_total-tracer_th_tend-tracer_src ! if (use_blobs) tracer_error = tracer_error + dmasstracer_L2 flux_error = tracer_error/(dtime*Grd%tcellsurf + epsln) if (T_prog(n)%name =='temp') then write (stdoutunit,'(/a)') ' ----Single time step diagnostics for tracer '//trim(T_prog(n)%name)//'----' write (stdoutunit,'(a,es24.17,a)') ' Total heat in ocean at (tau) referenced to 0degC = ',& ttracer,' J' if (use_blobs) then write (stdoutunit,'(a,es24.17,a)') ' Total heat change in L system for (taup1-taum1) = ',& dmasstracer_L,' J' write (stdoutunit,'(a,es24.17,a)') ' Total heat change in E system for (taup1-taum1) = ',& dmasstracer_E,' J' endif write (stdoutunit,'(a,es24.17,a)') ' Total heat change of full system for (taup1-taum1) = ',& dmasstracer,' J' if (use_blobs) then write (stdoutunit,'(a,es24.17,a)') ' Total heat transfer from L to E system (tau) = ',& dmasstracer_L2,' J' write (stdoutunit,'(a,es24.17,a)') ' Total heat in ocean L system at (tau) ref to 0degC = ',& tracer_L,' J' write (stdoutunit,'(a,es24.17,a)') ' Total heat in ocean E system at (tau) ref to 0degC = ',& tracer_E,' J' endif write (stdoutunit,'(a,es24.17,a)') ' Total heat input to ocean referenced to 0degC = ',& tracer_input_total,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via surface heat fluxes = ',& tracer_input_stf,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via bottom heat fluxes = ',& tracer_input_btf,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via open boundaries = ',& tracer_input_otf,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via precip-evap+melt = ',& tracer_input_pme,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via river runoff = ',& tracer_input_runoff,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via calving land ice = ',& tracer_input_calving,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via frazil formation = ',& frazil_heat,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via sources in the source array = ',& tracer_src,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via sources in th_tendency, or errors = ',& tracer_th_tend,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via eta_t smooth = ',& tracer_input_eta_smooth,' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via pbot_t smooth = ',& tracer_input_pbot_smooth,' J' write (stdoutunit,'(a,es24.17,a)') ' d(T*rho*dV) = ',& dmasstracer,' J' write (stdoutunit,'(a,es24.17,a)') ' Tracer mismatch: cp*d(rho*dV*T)-input = ',& tracer_error,' J' write (stdoutunit,'(a,es24.17,a)') ' Mismatch converted to a surface flux = ',& flux_error,' W/m^2' elseif (T_prog(n)%name(1:3) =='age') then tracer_error = dmasstracer-tracer_input_total-dtts*sec_in_yr_r*(mass_taup1-mass_taup1_keq1) tracer_error = tracer_error*mass_taup1_r flux_error = sec_in_yr*tracer_error/(dtime*Grd%tcellsurf + epsln) write (stdoutunit,'(/a)') ' ----Single time step diagnostics for tracer '//trim(T_prog(n)%name)//'----' write (stdoutunit,'(a,es24.17,a)') ' Total age at (tau) [sum(rho*dV*age)/mass_ocean] = ',& ttracer/mass_total,' yr' if (use_blobs) then if (mass_L .ne. 0.0) then write (stdoutunit,'(a,es24.17,a)') ' Total age in ocean L system at (tau) = ',& tracer_L/mass_L,' yr' else write (stdoutunit,'(a,es24.17,a)') ' Total age in ocean L system at (tau) = ',& 0.0,' yr' endif write (stdoutunit,'(a,es24.17,a)') ' Total age in ocean E system at (tau) = ',& tracer_E/mass_E,' yr' endif write (stdoutunit,'(a,es24.17,a)') ' dtts*(mass_ocean-mass_k=1)*sec_in_yr_r = ',& dtts*sec_in_yr_r*(mass_taup1-mass_taup1_keq1),' kg*yr' write (stdoutunit,'(a,es24.17,a)') ' dtts*(mass_ocean-mass_k=1)*sec_in_yr_r/mass_taup1 = ',& mass_taup1_r*dtts*sec_in_yr_r*(mass_taup1-mass_taup1_keq1),' yr' write (stdoutunit,'(a,es24.17,a)') ' Total age input to ocean = ',& tracer_input_total*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via surface fluxes = ',& tracer_input_stf*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via bottom fluxes = ',& tracer_input_btf*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via open boundaries = ',& tracer_input_otf*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via precip-evap+melt = ',& tracer_input_pme*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via river runoff = ',& tracer_input_runoff*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via calving land ice = ',& tracer_input_calving*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via sources in source array = ',& tracer_src*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via sources in th_tendency, or errors = ',& tracer_th_tend*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via eta_t smoother = ',& tracer_input_eta_smooth*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via pbot_t smoother = ',& tracer_input_pbot_smooth*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' d(age*rho*vol)/mass_ocean = ',& dmasstracer*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age mismatch: [d(rho*dV*age) - dtts*(Mt-M1)]/Mt = ',& tracer_error,' yr' write (stdoutunit,'(a,es24.17,a)') ' Mismatch converted to a surface flux = ',& flux_error,' 1/m^2' else write (stdoutunit,'(/a)') '----Single time step diagnostics for tracer '//trim(T_prog(n)%name)//'----' write (stdoutunit,'(a,es24.17,a)') ' Total tracer in ocean at time (tau) = ',& ttracer,' kg' if (use_blobs) then write (stdoutunit,'(a,es24.17,a)') ' Total tracer change in L system for (taup1-taum1) = ',& dmasstracer_L,' kg' write (stdoutunit,'(a,es24.17,a)') ' Total tracer change in E system for (taup1-taum1) = ',& dmasstracer_E,' kg' endif write (stdoutunit,'(a,es24.17,a)') ' Total tracer change in system for (taup1-taum1) = ',& dmasstracer,' kg' if (use_blobs) then write (stdoutunit,'(a,es24.17,a)') ' Total tracer transfer from L to E system (tau) = ',& dmasstracer_L2,' kg' write (stdoutunit,'(a,es24.17,a)') ' Total tracer in ocean L system at (tau) = ',& tracer_L,' kg' write (stdoutunit,'(a,es24.17,a)') ' Total tracer in ocean E system at (tau) = ',& tracer_E,' kg' endif write (stdoutunit,'(a,es24.17,a)') ' Total tracer input to ocean = ',& tracer_input_total,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via surface fluxes = ',& tracer_input_stf,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via bottom fluxes = ',& tracer_input_btf,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via open boundaries = ',& tracer_input_otf,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via precip-evap+melt = ',& tracer_input_pme,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via river runoff = ',& tracer_input_runoff,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via calving land ice = ',& tracer_input_calving,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via sources in source array = ',& tracer_src,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via sources in th_tendency, or errors = ',& tracer_th_tend,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via eta_t smoother = ',& tracer_input_eta_smooth,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via pbot_t smoother = ',& tracer_input_pbot_smooth,' kg' if (T_prog(n)%name =='salt') then write (stdoutunit,'(a,es24.17,a)') ' .001*d(S*rho*dV) = ', & psu2ppt*dmasstracer,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer mismatch: 0.001*d(rho*dV*T)-input = ', & psu2ppt*tracer_error,' kg ' write (stdoutunit,'(a,es24.17,a)') ' Mismatch converted to a surface flux = ', & psu2ppt*flux_error,' kg/(m^2 sec) ' else write (stdoutunit,'(a,es24.17,a)') ' d(T*rho*dV) = ',dmasstracer,' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer mismatch: d(rho*dV*T)-input = ',tracer_error,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mismatch converted to a surface flux = ',flux_error,' kg/(m^2 sec) ' endif endif enddo ! end n-loop call mpp_clock_end(id_tracer_change) end subroutine tracer_change ! NAME="tracer_change" !####################################################################### ! ! ! ! Compute integrated tracer in model. ! ! function total_tracer (Tracer, Thickness, index) type(ocean_prog_tracer_type), intent(in) :: Tracer type(ocean_thickness_type), intent(in) :: Thickness integer, intent(in) :: index real :: total_tracer real, dimension(isd:ied,jsd:jed) :: tracer_k integer :: k if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (total_tracer): module needs initialization ') endif total_tracer = 0.0 tracer_k(:,:) = 0.0 if (use_blobs) then do k=1,nk tracer_k(:,:) = tracer_k(:,:) + Tracer%conversion*Grd%tmask(:,:,k) & *( Grd%dat(:,:)*Thickness%rho_dzt(:,:,k,index)*Tracer%field(:,:,k,index) & + Tracer%sum_blob(:,:,k,index)) enddo else do k=1,nk tracer_k(:,:) = tracer_k(:,:) + Tracer%conversion*Grd%tmask(:,:,k)*Grd%dat(:,:) & *Thickness%rho_dzt(:,:,k,index)*Tracer%field(:,:,k,index) enddo endif if(have_obc) tracer_k(:,:) = tracer_k(:,:)*Grd%obc_tmask(:,:) total_tracer = mpp_global_sum(Dom%domain2d,tracer_k(:,:), global_sum_flag) end function total_tracer ! NAME="total_tracer" !####################################################################### ! ! ! ! Compute integrated tracer on a k-level. ! ! function klevel_total_tracer(Tracer, Thickness, tindex, kindex) type(ocean_prog_tracer_type), intent(in) :: Tracer type(ocean_thickness_type), intent(in) :: Thickness integer, intent(in) :: tindex integer, intent(in) :: kindex real :: klevel_total_tracer real, dimension(isd:ied,jsd:jed) :: tracer_k if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (total_tracer): module needs initialization ') endif klevel_total_tracer=0.0 if (use_blobs) then tracer_k(:,:) = Tracer%conversion*Grd%tmask(:,:,kindex) & *( Grd%dat(:,:)*Thickness%rho_dzt(:,:,kindex,tindex)*Tracer%field(:,:,kindex,tindex) & + Tracer %sum_blob(:,:,kindex,tindex) ) else tracer_k(:,:) = Tracer%conversion*Grd%tmask(:,:,kindex)*Grd%dat(:,:) & *Thickness%rho_dzt(:,:,kindex,tindex)*Tracer%field(:,:,kindex,tindex) endif if(have_obc) tracer_k(:,:) = tracer_k(:,:)*Grd%obc_tmask(:,:) klevel_total_tracer = mpp_global_sum(Dom%domain2d,tracer_k(:,:), global_sum_flag) end function klevel_total_tracer ! NAME="klevel_total_tracer" !####################################################################### ! ! ! ! Compute total ocean tracer cell mass. For Boussinesq fluid, ! mass is determined using rho0 for density. ! ! function total_mass (Thickness, index) type(ocean_thickness_type), intent(in) :: Thickness integer, intent(in) :: index real, dimension(isd:ied,jsd:jed) :: mass_t real :: total_mass integer :: k if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (total_mass): module needs initialization ') endif total_mass = 0.0 mass_t(:,:) = 0.0 if (use_blobs) then do k=1,nk mass_t(:,:) = mass_t(:,:) + Thickness%rho_dztT(:,:,k,index)*Grd%dat(:,:)*Grd%tmask(:,:,k) enddo else do k=1,nk mass_t(:,:) = mass_t(:,:) + Thickness%rho_dzt(:,:,k,index)*Grd%dat(:,:)*Grd%tmask(:,:,k) enddo endif if(have_obc) mass_t(:,:) = mass_t(:,:)*Grd%obc_tmask(:,:) total_mass = mpp_global_sum(Dom%domain2d,mass_t(:,:), global_sum_flag) end function total_mass ! NAME="total_mass" !####################################################################### ! ! ! ! Compute total ocean tracer cell volume. ! ! function total_volume (Thickness, index) type(ocean_thickness_type), intent(in) :: Thickness integer, intent(in) :: index real, dimension(isd:ied,jsd:jed) :: volume_t real :: total_volume integer :: k if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (total_volume): module needs initialization ') endif total_volume = 0.0 volume_t(:,:) = 0.0 if (use_blobs) then do k=1,nk volume_t(:,:) = volume_t(:,:) + Thickness%dztT(:,:,k,index)*Grd%dat(:,:)*Grd%tmask(:,:,k) enddo else do k=1,nk volume_t(:,:) = volume_t(:,:) + Thickness%dzt(:,:,k)*Grd%dat(:,:)*Grd%tmask(:,:,k) enddo endif if(have_obc) volume_t(:,:) = volume_t(:,:)*Grd%obc_tmask(:,:) total_volume = mpp_global_sum(Dom%domain2d, volume_t(:,:), global_sum_flag) end function total_volume ! NAME="total_volume" !####################################################################### ! ! ! ! Compute ocean tracer cell mass in a k-level. For Boussinesq fluid, ! mass is determined using rho0 for density. ! ! function klevel_total_mass (Thickness, tindex, kindex) type(ocean_thickness_type), intent(in) :: Thickness integer, intent(in) :: tindex integer, intent(in) :: kindex real, dimension(isd:ied,jsd:jed) :: mass_t real :: klevel_total_mass if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (klevel_total_mass): module needs initialization ') endif klevel_total_mass = 0.0 if (use_blobs) then mass_t(:,:) = Thickness%rho_dztT(:,:,kindex,tindex)*Grd%dat(:,:)*Grd%tmask(:,:,kindex) else mass_t(:,:) = Thickness%rho_dzt(:,:,kindex,tindex)*Grd%dat(:,:)*Grd%tmask(:,:,kindex) endif if(have_obc) mass_t(:,:) = mass_t(:,:)*Grd%obc_tmask(:,:) klevel_total_mass = mpp_global_sum(Dom%domain2d,mass_t(:,:), global_sum_flag) end function klevel_total_mass ! NAME="klevel_total_mass" !####################################################################### ! ! ! ! Compute some integrated tracer diagnostics. ! ! subroutine tracer_integrals (Time, Thickness, Tracer) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Tracer real, dimension(isd:ied,jsd:jed) :: mass_t real, dimension(isd:ied,jsd:jed,3) :: T_field real :: tbar, tvar, tchg, ttot real :: mass_total integer :: k, tau, taum1, taup1 integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (tracer_integrals): module needs initialization ') endif call mpp_clock_begin(id_tracer_integrals) taum1 = Time%taum1 tau = Time%tau taup1 = Time%taup1 write(stdoutunit,*) ' ' call write_timestamp(Time%model_time) write(stdoutunit,'(/a,a)') 'Tracer integrals for ' , trim(Tracer%name) tbar = 0.0; tvar = 0.0; tchg = 0.0; ttot = 0.0 mass_total = total_mass(Thickness, tau) wrk1_2d(:,:) = 0.0 wrk2_2d(:,:) = 0.0 wrk3_2d(:,:) = 0.0 if (use_blobs) then do k=1,nk mass_t(:,:) = Thickness%rho_dzt(:,:,k,tau)*Grd%dat(:,:)*Grd%tmask(:,:,k) if(have_obc) mass_t(:,:) = mass_t(:,:) * Grd%obc_tmask(:,:) wrk1_2d(:,:) = wrk1_2d(:,:) + Tracer%field(:,:,k,tau)*mass_t(:,:) & + Tracer%sum_blob(:,:,k,tau)*Grd%tmask(:,:,k) T_field(:,:,tau) = ( Grd%datr(:,:)*Tracer%sum_blob(:,:,k,tau) & +Tracer%field(:,:,k,tau)*Thickness%rho_dzt(:,:,k,tau) )& /(Thickness%rho_dztT(:,:,k,tau) + epsln) wrk2_2d(:,:) = wrk2_2d(:,:) + T_field(:,:,tau)*T_field(:,:,tau)& *Thickness%rho_dztT(:,:,k,tau)*Grd%dat(:,:) T_field(:,:,taup1) = ( Grd%datr(:,:)*Tracer%sum_blob(:,:,k,taup1) & +Tracer%field(:,:,k,taup1)*Thickness%rho_dzt(:,:,k,taup1) )& /(Thickness%rho_dztT(:,:,k,taup1) + epsln) T_field(:,:,taum1) = ( Grd%datr(:,:)*Tracer%sum_blob(:,:,k,taum1) & +Tracer%field(:,:,k,taum1)*Thickness%rho_dzt(:,:,k,taum1) )& /(Thickness%rho_dztT(:,:,k,taum1) + epsln) wrk3_2d(:,:) = wrk3_2d(:,:) + abs(T_field(:,:,taup1)-T_field(:,:,taum1))*dtimer enddo else do k=1,nk mass_t(:,:) = Thickness%rho_dzt(:,:,k,tau)*Grd%dat(:,:)*Grd%tmask(:,:,k) if(have_obc) mass_t(:,:) = mass_t(:,:) * Grd%obc_tmask(:,:) wrk1_2d(:,:) = wrk1_2d(:,:) + Tracer%field(:,:,k,tau)*mass_t(:,:) wrk2_2d(:,:) = wrk2_2d(:,:) + Tracer%field(:,:,k,tau)*Tracer%field(:,:,k,tau)*mass_t(:,:) wrk3_2d(:,:) = wrk3_2d(:,:) + abs(Tracer%field(:,:,k,taup1)-Tracer%field(:,:,k,taum1))*dtimer enddo endif tbar = mpp_global_sum(Dom%domain2d, wrk1_2d(:,:), global_sum_flag) tvar = mpp_global_sum(Dom%domain2d, wrk2_2d(:,:), global_sum_flag) tchg = mpp_global_sum(Dom%domain2d, wrk3_2d(:,:), global_sum_flag) ttot = total_tracer(Tracer,Thickness,tau) tbar = tbar/mass_total tvar = tvar/mass_total - tbar**2 write(stdoutunit, '(/a,a,a,es24.17)') 'Average ',trim(Tracer%name),' = ', tbar write(stdoutunit, '(a,a,a,es24.17)') 'Variance ',trim(Tracer%name),' = ', tvar write(stdoutunit, '(a,a,a,es24.17)') '|dT/dt| ',trim(Tracer%name),' = ', tchg write(stdoutunit, '(a,a,a,es24.17/)') 'Total ',trim(Tracer%name),' = ', ttot call mpp_clock_end(id_tracer_integrals) end subroutine tracer_integrals ! NAME="tracer_integrals" !####################################################################### ! ! ! ! Check to be sure ocean tracer is zero over land ! ! subroutine tracer_land_cell_check (Time, T_prog) type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), intent(in) :: T_prog(:) integer :: i, j, k, n, num, taup1 integer :: stdoutunit stdoutunit=stdout() call mpp_clock_begin(id_tracer_land_cell_check) if(size(T_prog,1) /= num_prog_tracers) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (tracer_land_cell_check): passing T_prog with wrong dimensions') endif if (.not. module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag_mod (tracer_land_cell_check): module needs initialization ') endif taup1 = Time%taup1 write (stdoutunit,'(1x,/a/)')'Locations (if any) where land cell tracer is non-zero...' num = 0 do n=1,num_prog_tracers do k=1,nk do j=jsc,jec do i=isc,iec if (Grd%tmask(i,j,k) == 0 .and. (T_prog(n)%field(i,j,k,taup1) /= 0.0) ) then num = num + 1 if (num < 10) then write(unit,9100) & i+Dom%ioff,j+Dom%joff,k,Grd%xt(i,j),Grd%yt(i,j),Grd%zt(k),n,T_prog(n)%field(i,j,k,taup1) endif endif enddo enddo enddo enddo if (num > 0) call mpp_error(FATAL) 9100 format(/' " =>Error: Land cell at (i,j,k) = ','(',i4,',',i4,',',i4,'),',& ' (lon,lat,dpt) = (',f8.3,',',f8.3,',',f12.3,'m) has tracer(n=',i2,') =',e10.3) call mpp_clock_end(id_tracer_land_cell_check) end subroutine tracer_land_cell_check ! NAME="tracer_land_cell_check" !####################################################################### ! ! ! ! Compute change in mass over many time steps, and compare to the ! input of mass through surface to check for mass conservation. ! !============================================================ ! ! threelevel scheme ! ! Here is the logic for the accumulation of the fluxes and ! comparisons between mass/volumes at the start and the end. ! ! Consider accumulation over four leap-frog time steps. ! Ignore time filtering. ! ! mass(2) = mass(0) + 2dt*F(1) taup1=2, taum1=0, tau=1 ! ! mass(3) = mass(1) + 2dt*F(2) taup1=3, taum1=1, tau=2 ! ! mass(4) = mass(2) + 2dt*F(3) taup1=4, taum1=2, tau=3 ! ! mass(5) = mass(3) + 2dt*F(4) taup1=5, taum1=3, tau=4 ! ! Hence, ! ! [mass(4) + mass(5)] = [mass(0) + mass(1)] + 2dt*[F(1)+F(2)+F(3)+F(4)] ! ! For this example, we have ! ! itts_mass=1 through itte_mass=4 for accumulating fluxes ! ! itt=itts_mass=1=tau we use taum1=0 and tau=1 to get starting mass ! ! itt=itte_mass=4=tau we use tau=4 and taup1=5 to get the final mass ! !============================================================ ! ! twolevel scheme ! ! Here is the logic for the accumulation of the fluxes and ! comparisons between mass/volumes at the start and the end. ! ! Consider accumulation over four time steps. ! ! mass(3/2) = mass(1/2) + dt*F(1) taup1=3/2, taum1=1/2, tau=1 ! ! mass(5/2) = mass(3/2) + dt*F(2) taup1=5/2, taum1=3/2, tau=2 ! ! mass(7/2) = mass(5/2) + dt*F(3) taup1=7/2, taum1=5/2, tau=3 ! ! mass(9/2) = mass(7/2) + dt*F(4) taup1=9/2, taum1=7/2, tau=4 ! ! Hence, ! ! mass(9/2) = mass(1/2) + dt*[F(1)+F(2)+F(3)+F(4)] ! ! For this example, we have ! ! itts_mass=1 through itte_mass=4 for accumulating fluxes ! ! itt=itts_mass=1=tau we use taum1=1/2 to get starting mass ! ! itt=itte_mass=4=tau we use taup1=9/2 to get the final mass ! ! ! subroutine mass_conservation (Time, Thickness, Ext_mode, pme, runoff, calving) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_external_mode_type), intent(in) :: Ext_mode real, dimension(isd:,jsd:), intent(in) :: pme real, dimension(isd:,jsd:), intent(in) :: runoff real, dimension(isd:,jsd:), intent(in) :: calving real :: total_time real :: mass_error, mass_error_rate real :: mass_input, mass_eta_smooth, mass_pbot_smooth, mass_chg real :: pme_input, runoff_input, calving_input, source_input integer :: i, j integer :: itt, taum1, tau, taup1 logical, save :: first=.true. integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (mass_conservation): module needs initialization ') endif itt = Time%itt taum1 = Time%taum1 tau = Time%tau taup1 = Time%taup1 if (first) then first = .false. allocate (pme_flux(isd:ied,jsd:jed)) allocate (runoff_flux(isd:ied,jsd:jed)) allocate (source_flux(isd:ied,jsd:jed)) allocate (calving_flux(isd:ied,jsd:jed)) allocate (eta_smooth(isd:ied,jsd:jed)) allocate (pbot_smooth(isd:ied,jsd:jed)) pme_flux = 0.0 runoff_flux = 0.0 source_flux = 0.0 calving_flux = 0.0 eta_smooth = 0.0 pbot_smooth = 0.0 mass_start = 0.0 mass_final = 0.0 itts_mass = itt+2 itte_mass = itts_mass-2 + nint(tracer_conserve_days*86400.0/dtts) if(itts_mass >= itte_mass) then write(stdoutunit,*) & '==>Warning: increase tracer_conserve_days to properly compute mass and tracer conservation.' endif endif if(tendency==THREE_LEVEL) then if (itt == itts_mass) then mass_start = 0.5*(total_mass(Thickness, taum1)+total_mass(Thickness, tau)) endif if (itt == itte_mass) then mass_final = 0.5*(total_mass(Thickness, tau)+total_mass(Thickness, taup1)) endif elseif(tendency==TWO_LEVEL) then if (itt == itts_mass) then mass_start = total_mass(Thickness, taum1) endif if (itt == itte_mass) then mass_final = total_mass(Thickness, taup1) endif endif if(itt >= itts_mass .and. itt <= itte_mass) then do j=jsc,jec do i=isc,iec pme_flux(i,j) = pme_flux(i,j) + dteta*Grd%dat(i,j)*pme(i,j) runoff_flux(i,j) = runoff_flux(i,j) + dteta*Grd%dat(i,j)*runoff(i,j) calving_flux(i,j) = calving_flux(i,j)+ dteta*Grd%dat(i,j)*calving(i,j) source_flux(i,j) = source_flux(i,j) + dteta*Grd%dat(i,j)*Ext_mode%source(i,j) eta_smooth(i,j) = eta_smooth(i,j) + dteta*Grd%tmask(i,j,1)*Grd%dat(i,j)*Ext_mode%eta_smooth(i,j) pbot_smooth(i,j) = pbot_smooth(i,j) + dteta*Grd%tmask(i,j,1)*Grd%dat(i,j)*Ext_mode%pbot_smooth(i,j) enddo enddo endif if (itt==itte_mass) then if(have_obc) then pme_flux(:,:) = pme_flux(:,:)*Grd%obc_tmask(:,:) runoff_flux(:,:) = runoff_flux(:,:)*Grd%obc_tmask(:,:) calving_flux(:,:) = calving_flux(:,:)*Grd%obc_tmask(:,:) source_flux(:,:) = source_flux(:,:)*Grd%obc_tmask(:,:) eta_smooth(:,:) = eta_smooth(:,:)*Grd%obc_tmask(:,:) pbot_smooth(:,:) = pbot_smooth(:,:)*Grd%obc_tmask(:,:) eta_smooth(:,:) = eta_smooth(:,:)*Grd%obc_tmask(:,:) pbot_smooth(:,:) = pbot_smooth(:,:)*Grd%obc_tmask(:,:) endif pme_input = mpp_global_sum(Dom%domain2d,pme_flux(:,:),global_sum_flag) runoff_input = mpp_global_sum(Dom%domain2d,runoff_flux(:,:),global_sum_flag) calving_input = mpp_global_sum(Dom%domain2d,calving_flux(:,:),global_sum_flag) source_input = mpp_global_sum(Dom%domain2d,source_flux(:,:),global_sum_flag) mass_eta_smooth = mpp_global_sum(Dom%domain2d,eta_smooth(:,:), global_sum_flag) mass_pbot_smooth = mpp_global_sum(Dom%domain2d,pbot_smooth(:,:),global_sum_flag) mass_input = pme_input + runoff_input + calving_input + source_input + mass_eta_smooth + mass_pbot_smooth total_time = (itte_mass-itts_mass+1)*dtts+epsln mass_chg = mass_final - mass_start mass_error = mass_chg-mass_input mass_error_rate = mass_error/total_time write (stdoutunit,'(/50x,a/)') & ' Measures of global integrated mass conservation over multiple time steps' write (stdoutunit,'(a,i10,a,es24.17,a)') & ' Ocean mass at timestep ',itts_mass,' = ',mass_start,' kg.' write (stdoutunit,'(a,i10,a,es24.17,a)') & ' Ocean mass at timestep ',itte_mass,' = ',mass_final,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Change in ocean mass over time interval = ',mass_chg,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Mass input via all contributions over time interval = ',mass_input,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Error in mass content change = ',mass_error,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Error in rate of mass content change = ',mass_error_rate,' kg/s.' write (stdoutunit,'(a,es24.17,a)') & ' Mass input via P-E fluxes over time interval = ',pme_input,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Mass input via runoff fluxes over time interval = ',runoff_input,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Mass input via ice calving fluxes over time interval = ',calving_input,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Mass input via source fluxes over time interval = ',source_input,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Mass input via eta_t smoother over time interval = ',mass_eta_smooth,' kg.' write (stdoutunit,'(a,es24.17,a)') & ' Mass input via pbot_t smoother over time interval = ',mass_pbot_smooth,' kg.' write (stdoutunit,'(/)') endif end subroutine mass_conservation ! NAME="mass_conservation" !####################################################################### ! ! ! ! Compute change in global integrated tracer over many time steps, ! and compare to the input of tracer through the boundaries to ! check for total tracer conservation. ! ! Accumulate fluxes as in the mass_conservation diagnostic. ! ! ! subroutine tracer_conservation (Time, Thickness, T_prog, T_diag, pme, runoff, calving) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: T_prog(:) type(ocean_diag_tracer_type), intent(in) :: T_diag(:) real, dimension(isd:,jsd:), intent(in) :: pme real, dimension(isd:,jsd:), intent(in) :: runoff real, dimension(isd:,jsd:), intent(in) :: calving real, dimension(isd:ied,jsd:jed) :: tmp real, dimension(isd:ied,jsd:jed) :: dt_con_area real, dimension(isd:ied,jsd:jed) :: runoff_mod real :: temporary real :: tracer_chg real :: tracer_error real :: tracer_error_rate real :: tracer_tend_input real :: tracer_calving_input real :: tracer_runoff_input real :: tracer_pme_input real :: tracer_stf_input real :: tracer_btf_input real :: tracer_otf_input real :: tracer_frazil_input real :: tracer_total_input real :: tracer_source_input real :: tracer_eta_smooth_input real :: tracer_pbot_smooth_input real :: total_time integer :: i, j, k, n integer :: itt integer :: tau, taum1, taup1 logical, save :: first =.true. integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (tracer_conservation): module needs initialization ') endif if(size(T_prog,1) /= num_prog_tracers) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (tracer_conservation): passing T_prog with wrong dimensions') endif tau = Time%tau taum1 = Time%taum1 taup1 = Time%taup1 itt = Time%itt if (first) then first = .false. allocate (tracer_tend(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_runoff(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_calving(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_pme(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_otf(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_stf(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_btf(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_frazil(isd:ied,jsd:jed)) allocate (tracer_source(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_eta_smooth(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_pbot_smooth(isd:ied,jsd:jed,num_prog_tracers)) allocate (tracer_start(num_prog_tracers)) allocate (tracer_final(num_prog_tracers)) tracer_source(:,:,:) = 0.0 tracer_tend(:,:,:) = 0.0 tracer_runoff(:,:,:) = 0.0 tracer_calving(:,:,:) = 0.0 tracer_pme(:,:,:) = 0.0 tracer_otf(:,:,:) = 0.0 tracer_stf(:,:,:) = 0.0 tracer_btf(:,:,:) = 0.0 tracer_frazil(:,:) = 0.0 tracer_source(:,:,:) = 0.0 tracer_eta_smooth(:,:,:) = 0.0 tracer_pbot_smooth(:,:,:) = 0.0 tracer_start(:) = 0.0 tracer_final(:) = 0.0 itts_tracer = itt+4 itte_tracer = itts_tracer-4 + nint(tracer_conserve_days*86400.0/dtts) endif if(tendency==THREE_LEVEL) then if (itt == itts_tracer) then do n=1,num_prog_tracers tracer_start(n) = 0.5*( total_tracer(T_prog(n),Thickness,taum1) & + total_tracer(T_prog(n),Thickness,tau)) enddo endif if (itt == itte_tracer) then do n=1,num_prog_tracers tracer_final(n) = 0.5*( total_tracer(T_prog(n),Thickness,tau) & + total_tracer(T_prog(n),Thickness,taup1)) enddo endif elseif(tendency==TWO_LEVEL) then if (itt == itts_tracer) then do n=1,num_prog_tracers tracer_start(n) = total_tracer(T_prog(n),Thickness,taum1) enddo endif if (itt == itte_tracer) then do n=1,num_prog_tracers tracer_final(n) = total_tracer(T_prog(n),Thickness,taup1) enddo endif endif if(itt >= itts_tracer .and. itt <= itte_tracer) then ! runoff can be negative, as for the "geoengineering" trick of siphoning ! Arabian Sea water to the Caspian Sea, to keep the Caspian from evaporating ! away during a coupled climate simulation. As negative runoff is ignored ! in the rivermix module (rightly so, as we do not wish to insert tracer ! from negative runoff into the ocean), we also wish to ignore negative ! runoff for the diagnostics here, in order to reflect what the prognostic ! equations are doing. runoff_mod(:,:) = 0.0 do j=jsc,jec do i=isc,iec if(runoff(i,j) > 0) then runoff_mod(i,j) = runoff(i,j) endif enddo enddo do n=1,num_prog_tracers ! global integral of th_tendency leaves just the sources ! (and other stuff to be removed below) ! define an intermediate array do j=jsc,jec do i=isc,iec dt_con_area(i,j) = dtts*T_prog(n)%conversion*Grd%dat(i,j) enddo enddo do k=1,nk do j=jsc,jec do i=isc,iec tracer_tend(i,j,n) = tracer_tend(i,j,n) + & dt_con_area(i,j)*Grd%tmask(i,j,k)*T_prog(n)%th_tendency(i,j,k) tracer_source(i,j,n) = tracer_source(i,j,n) + & dt_con_area(i,j)*Grd%tmask(i,j,k)*Thickness%rho_dzt(i,j,k,taup1)*T_prog(n)%source(i,j,k) enddo enddo enddo ! btf < 0 means heat is added to ocean; hence the minus sign do j=jsc,jec do i=isc,iec temporary = dt_con_area(i,j)*Grd%tmask(i,j,1) tracer_eta_smooth(i,j,n) = tracer_eta_smooth(i,j,n) + temporary*T_prog(n)%eta_smooth(i,j) tracer_pbot_smooth(i,j,n) = tracer_pbot_smooth(i,j,n) + temporary*T_prog(n)%pbot_smooth(i,j) tracer_runoff(i,j,n) = tracer_runoff(i,j,n) + temporary*T_prog(n)%trunoff(i,j)*runoff_mod(i,j) tracer_calving(i,j,n) = tracer_calving(i,j,n) + temporary*T_prog(n)%tcalving(i,j)*calving(i,j) tracer_pme(i,j,n) = tracer_pme(i,j,n) + temporary*T_prog(n)%tpme(i,j)*pme(i,j) tracer_stf(i,j,n) = tracer_stf(i,j,n) + temporary*T_prog(n)%stf(i,j) tracer_btf(i,j,n) = tracer_btf(i,j,n) - temporary*T_prog(n)%btf(i,j) enddo enddo if(n==index_temp .and. index_frazil > 0) then do k=1,nk do j=jsc,jec do i=isc,iec tracer_frazil(i,j) = tracer_frazil(i,j) & + T_diag(index_frazil)%field(i,j,k)*Grd%dat(i,j)*Grd%tmask(i,j,k) enddo enddo enddo endif if (have_obc) then call ocean_obc_tracer_flux(T_prog(n), tmp) tracer_otf(:,:,n) = (tracer_otf(:,:,n) + tmp(:,:))*Grd%obc_tmask(:,:) tracer_stf(:,:,n) = tracer_stf(:,:,n)*Grd%obc_tmask(:,:) tracer_btf(:,:,n) = tracer_btf(:,:,n)*Grd%obc_tmask(:,:) tracer_runoff(:,:,n) = tracer_runoff(:,:,n)*Grd%obc_tmask(:,:) tracer_calving(:,:,n) = tracer_calving(:,:,n)*Grd%obc_tmask(:,:) tracer_pme(:,:,n) = tracer_pme(:,:,n)*Grd%obc_tmask(:,:) tracer_tend(:,:,n) = tracer_tend(:,:,n)*Grd%obc_tmask(:,:) tracer_source(:,:,n) = tracer_source(:,:,n)*Grd%obc_tmask(:,:) tracer_eta_smooth(:,:,n) = tracer_eta_smooth(:,:,n)*Grd%obc_tmask(:,:) tracer_pbot_smooth(:,:,n) = tracer_pbot_smooth(:,:,n)*Grd%obc_tmask(:,:) if (n==index_temp) then tracer_frazil(:,:) = tracer_frazil(:,:)*Grd%obc_tmask(:,:) endif endif enddo endif if (itt==itte_tracer) then write (stdoutunit,'(/50x,a/)') & ' Measures of global integrated tracer conservation over multiple time steps' total_time = (itte_mass-itts_mass+1)*dtts+epsln do n=1,num_prog_tracers tracer_frazil_input = 0.0 if (n==index_temp) then tracer_frazil_input = mpp_global_sum(Dom%domain2d,tracer_frazil(:,:), global_sum_flag) endif tracer_runoff_input = mpp_global_sum(Dom%domain2d,tracer_runoff(:,:,n), global_sum_flag) tracer_calving_input = mpp_global_sum(Dom%domain2d,tracer_calving(:,:,n), global_sum_flag) tracer_pme_input = mpp_global_sum(Dom%domain2d,tracer_pme(:,:,n), global_sum_flag) tracer_otf_input = mpp_global_sum(Dom%domain2d,tracer_otf(:,:,n), global_sum_flag) tracer_stf_input = mpp_global_sum(Dom%domain2d,tracer_stf(:,:,n), global_sum_flag) tracer_btf_input = mpp_global_sum(Dom%domain2d,tracer_btf(:,:,n), global_sum_flag) tracer_tend_input = mpp_global_sum(Dom%domain2d,tracer_tend(:,:,n), global_sum_flag) tracer_source_input = mpp_global_sum(Dom%domain2d,tracer_source(:,:,n), global_sum_flag) tracer_eta_smooth_input = mpp_global_sum(Dom%domain2d,tracer_eta_smooth(:,:,n), global_sum_flag) tracer_pbot_smooth_input = mpp_global_sum(Dom%domain2d,tracer_pbot_smooth(:,:,n), global_sum_flag) tracer_total_input = tracer_stf_input + tracer_btf_input + tracer_otf_input & + tracer_runoff_input + tracer_calving_input + tracer_pme_input & + tracer_frazil_input + tracer_eta_smooth_input + tracer_pbot_smooth_input ! runoff, calving, pme, and smooth are added to T_prog(n)%th_tendency inside ! ocean_rivermix_mod (runoff, calving) and ocean_tracer_mod (pme, and smooth). ! To avoid double counting their contribution, we remove these terms from ! tracer_tend_input. As a result, tracer_tend_input will only include ! contributions from the global integration of internal flux convergences. ! When integrated globally, these convergences should sum to zero, and such is ! an important check that the contributions to tracer updates are coded properly. tracer_tend_input = tracer_tend_input - tracer_runoff_input - tracer_calving_input & -tracer_pme_input - tracer_eta_smooth_input - tracer_pbot_smooth_input ! stf and btf added to th_tendency when use explicit vertical diffusion (aidif=0) if (aidif==0.0) tracer_tend_input = tracer_tend_input - tracer_stf_input - tracer_btf_input tracer_chg = tracer_final(n)-tracer_start(n) tracer_error = tracer_chg - tracer_total_input - tracer_tend_input - tracer_source_input tracer_error_rate = tracer_error/(total_time*Grd%tcellsurf) if (n==index_temp) then write (stdoutunit,'(a,i10,a,es24.17,a)') ' Ocean heat content at timestep ',itts_tracer,' = ',tracer_start(n),' J.' write (stdoutunit,'(a,i10,a,es24.17,a)') ' Ocean heat content at timestep ',itte_tracer,' = ',tracer_final(n),' J.' write (stdoutunit,'(a,es24.17,a)') ' Change in ocean heat content = ',tracer_chg,' J.' write (stdoutunit,'(a,es24.17,a)') ' Total Heat input = ',tracer_total_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Error in heat content change = ',tracer_error,' J.' write (stdoutunit,'(a,es24.17,a)') ' Error in rate of heat content change = ',tracer_error_rate,' W/m^2.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by tendency terms = ',tracer_tend_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by surface fluxes = ',tracer_stf_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by bottom fluxes = ',tracer_btf_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by runoff and calving = ',tracer_runoff_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by precip-evap+calving = ',tracer_pme_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by open boundaries = ',tracer_otf_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by frazil formation = ',tracer_frazil_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by eta_t smoother = ',tracer_eta_smooth_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by pbot_t smoother = ',tracer_pbot_smooth_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by sources = ',tracer_source_input,' J.' write (stdoutunit,'(/)') elseif(T_prog(n)%name(1:3)=='age') then write (stdoutunit,'(a,i10,a,es24.17,a)') ' Mass weighted '//trim(T_prog(n)%name)// & ' at timestep ',itts_tracer,' = ',tracer_start(n),' yr*kg.' write (stdoutunit,'(a,i10,a,es24.17,a)') ' Mass weighted '//trim(T_prog(n)%name)// & ' at timestep ',itte_tracer,' = ',tracer_final(n),' yr*kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' total tracer input = ',tracer_total_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' change in ocean tracer = ',tracer_chg,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' error in tracer change = ',tracer_error,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' error in rate of tracer change = ',sec_in_yr*tracer_error_rate,' kg/(m^2).' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by tendency terms = ',tracer_tend_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by surface fluxes = ',tracer_stf_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by bottom fluxes = ',tracer_btf_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by runoff and calving = ',tracer_runoff_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by precip-evap+calving = ',tracer_pme_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by open boundaries = ',tracer_otf_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by eta_t smoother = ',tracer_eta_smooth_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by pbot_t smoother = ',tracer_pbot_smooth_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by sources = ',tracer_source_input,' kg*yr.' write (stdoutunit,'(/)') else write (stdoutunit,'(a,i10,a,es24.17,a)') ' Total '//trim(T_prog(n)%name)// & ' mass at timestep ',itts_tracer,' = ',tracer_start(n),' kg' write (stdoutunit,'(a,i10,a,es24.17,a)') ' Total '//trim(T_prog(n)%name)// & ' mass at timestep ',itte_tracer,' = ',tracer_final(n),' kg' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' change in ocean tracer mass = ',tracer_chg,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' error in mass change = ',tracer_error,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' error in rate of mass change = ',tracer_error_rate,' kg/(m^2 sec).' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by tendency terms = ',tracer_tend_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by surface fluxes = ',tracer_stf_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by bottom fluxes = ',tracer_btf_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by runoff and calving = ',tracer_runoff_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by precip-evap+calving = ',tracer_pme_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by open boundaries = ',tracer_otf_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by eta_t smoother = ',tracer_eta_smooth_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by pbot_t smoother = ',tracer_pbot_smooth_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by sources = ',tracer_source_input,' kg.' write (stdoutunit,'(/)') endif enddo endif end subroutine tracer_conservation ! NAME="tracer_conservation" !####################################################################### ! ! ! ! Routine to diagnose the amount of mixing between classes of a ! particular tracer. Temperature is used as default. ! Method follows that used in the paper ! ! Spurious diapycnal mixing associated with advection in a ! z-coordinate ocean model, 2000: S.M. Griffies, R.C. ! Pacanowski, and R.W. Hallberg. Monthly Weather Review, vol 128, 538--564. ! ! This diagnostic is most useful when computing the levels of ! effective dia-tracer mixing occuring in a model run with ! zero buoyancy forcing at the boundaries. ! ! Algorithm notes: ! ! -assumes flat ocean bottom--non-flat bottoms loose the precise relation ! between sorted depth and true ocean depth. This is a minor inconvenience. ! The code is actually written so that the horizontal area of each ! layer can be different. This will allow for this diagnostic to be ! used, say, for simple topography, such as bowl or bump. ! ! -assumes Boussinesq fluid so that consider volume instead of mass of a cell. ! Also his means that dzt = rho0r*rho_dzt ! ! -assumes area integrated eta_t is zero, so domain volume is static. ! This is the case when there are no water boundary fluxes. ! ! -Results are meaningful only when dxt*dyt*dst of each grid cell ! is the same. This is best realized with a beta-plane or f-plane ! geometry, and with zstar vertical coordinate. ! ! My best understanding of this limitation is related to systematic ! biases in roundoff errors that result when the grid cells have ! varying volumes. ! ! If choose to use geopotential vertical coordinate, it is best ! to set linear_free_surface=.true. in ocean_thickness_nml, ! so that Thickness%rho_dzt = rho0%Grd%dzt. The sorting model of ! mixing has not been generalized to evolving layer thicknesses ! with geopotential. ! ! With zstar, the dst is constant in time, and the sorting method ! will sort to a depth in zstar space rather than geopotential space. ! This is a trivial distinction in principle, but should help with ! some roundoff issues in practice. ! ! -assumes tendency=TWO_LEVEL, which is exploited here to save memory. ! ! Numerical roundoff is a real issue with this diagnostic. ! It is critical that full double precision be used ! to garner sensible results. ! ! -defines some global arrays, so requires large memory. ! This feature can be removed if parallel sort is ! implemented. So far, such has not been done. ! ! -Effective kappa is set to zero at top of top-most cell. ! It is then diagnosed as zero (or roundoff) at bottom ! of the column if there are zero boundary buoyancy fluxes. ! ! -when computing density, we do rho=-alpha*theta. ! We drop the rho0 factor in order to reduce roundoff. ! Likewise, we assume alpha=1.0 rather than alpha=alpha_linear_eos ! We use alpha=1.0 to improve precision. ! With alpha=1.0 and rho=-alpha*theta, the rho variable ! is then just minus theta. ! ! -minimum vertical density gradient rho_grad_min is necessary to ! avoid errors with truncation in the division by drho/dz when compute kappa. ! rho_grad_min corresponds roughly to the precision of the computation. ! Physically, with ! ! N^2 = -(g/rho0)(drho/dz) ! ! then rho_grad_min sets a minimum N^2 resolved. ! This corresponds to a frequency f=N/2pi. The typical ! period of inertial oscillations in the deep ocean is 6hrs ! (Pickard and Emery, page 55-56). In the upper ocean, it is ! 10-30 minutes, in pycnocline it is smaller still. ! So to cover the majority of the ocean's stratification, ! we will want to set rho_grad_min to something smaller than 9e-6. ! ! To bin the effective diffusivity, it is also useful to have a max ! vertical density gradient. ! !-versions: ! ! mom4p0 method assumed rigid lid, or zero surface height undulations. ! Fit the following equation to the model data ! \partial_{t}(rho_sort) = (F_{k}-F_{k-1})/dzw ! ! mom4p1 method fits the following equation to model data ! \partial_{t}(dzt_sort*theta_sort) = F_{k} - F_{k-1} ! Fits this equation assuming two-time level tendency ! ! revision: 05/2005 ! revision: 07/2007 ! Stephen.Griffies ! ! ! subroutine diagnose_kappa_sort(Time, Thickness, Theta) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Theta real, dimension(:,:,:), allocatable :: global_tmask ! for global mask real, dimension(:,:,:), allocatable :: global_dzt ! for global dzt real, dimension(:,:,:), allocatable :: global_tracer ! for global tracer real, dimension(:,:) , allocatable :: global_dat ! for global area real, dimension(:) , allocatable :: rho_sort ! density (kg/m^3) of sorted parcels real, dimension(:) , allocatable :: vol_sort ! volume (m^3) of sorted parcels logical, save :: first_enter_routine=.true. integer :: tau, taup1 integer :: i, j, k, m, n integer :: nsort, numzero real :: kappa_prev, tmp real :: difference integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (diagnose_kappa_sort): module needs initialization ') endif ! taum1=tau since tendency=TWO_LEVEL tau = Time%tau taup1 = Time%taup1 allocate (global_tmask(ni,nj,nk)) ; global_tmask=0.0 call mpp_global_field(Dom%domain2d, Grd%tmask, global_tmask) allocate (global_dat(ni,nj)) ; global_dat=0.0 call mpp_global_field(Dom%domain2d, Grd%dat, global_dat) if (first_enter_routine) then first_enter_routine = .false. nsortpts = Grd%wet_t_points thick_column= 0.0 mass_column = 0.0 ! for diagnosing kappa_sort allocate (nlayer(nk)) ; nlayer = 0 allocate (wetarea(nk)) ; wetarea = 0.0 allocate (normalize(nk)) ; normalize = 0.0 allocate (deriv_lay(nk)) ; deriv_lay = 0.0 allocate (flux_lay(nk)) ; flux_lay = 0.0 allocate (kappa_lay(nk)) ; kappa_lay = 0.0 allocate (rho_lay(nk,2)) ; rho_lay = 0.0 allocate (dzt_rho_lay(nk,2)) ; dzt_rho_lay = 0.0 allocate (dzt_lay(nk,2)) ; dzt_lay = 0.0 allocate (dzw_lay(nk,2)) ; dzw_lay = 0.0 allocate (mass_lay(nk,2)) ; mass_lay = 0.0 allocate (volume_lay(nk,2)) ; volume_lay = 0.0 do k=1,nk nlayer(k) = 0 wetarea(k)= 0.0 do j=1,nj do i=1,ni if(global_tmask(i,j,k) == 1.0) then nlayer(k) = nlayer(k) + 1 wetarea(k) = wetarea(k) + global_dat(i,j) endif enddo enddo if(nlayer(k) > 0) then normalize(k) = rho0r/wetarea(k) else normalize(k) = 0.0 endif enddo else do k=1,nk rho_lay(k,1) = rho_lay(k,2) dzt_rho_lay(k,1) = dzt_rho_lay(k,2) dzt_lay(k,1) = dzt_lay(k,2) dzw_lay(k,1) = dzw_lay(k,2) mass_lay(k,1) = mass_lay(k,2) volume_lay(k,1) = volume_lay(k,2) enddo endif allocate (global_dzt(ni,nj,nk)) ; global_dzt=0.0 allocate (global_tracer(ni,nj,nk)) ; global_tracer=0.0 allocate (vol_sort(nsortpts)) ; vol_sort=0.0 allocate (rho_sort(nsortpts)) ; rho_sort=0.0 ! taup1 time level ! remove rho0 from dzt at a later point call mpp_global_field(Dom%domain2d, Theta%field(:,:,:,taup1), global_tracer) call mpp_global_field(Dom%domain2d, Thickness%dst(:,:,:), global_dzt) nsort=0 do k=1,nk do j=1,nj do i=1,ni if(global_tmask(i,j,k) == 1.0) then nsort=nsort+1 vol_sort(nsort) = rho0*global_dzt(i,j,k)*global_dat(i,j) rho_sort(nsort) = -global_tracer(i,j,k) ! assume alpha=1.0 endif enddo enddo enddo deallocate(global_tracer) deallocate(global_dzt) deallocate(global_dat) if(debug_diagnose_mixingB) then write(stdoutunit,'(/a)') '==========Before sort============' do nsort=1,nsortpts write(stdoutunit,'(a,i7,a,e17.11,a,i7,a,e17.11)') & ' rho_sort(',nsort,') = ',rho_sort(nsort), & ' vol_sort(',nsort,') = ',vol_sort(nsort)*rho0r enddo endif ! reorder rho_sort and carry vol_sort along with it. ! after the sorting, we will have ! rho_sort(1) = smallest rho_sort value, and vol_sort(1) is its volume ! rho_sort(nsortpts) = largest rho_sort value, and vol_sort(nsortpts) is its volume call sort_shell_array(rho_sort,vol_sort) if(debug_diagnose_mixingC) then write(stdoutunit,'(/a)') '==========After sort============' do nsort=1,nsortpts write(stdoutunit,'(a,i7,a,e17.11,a,i7,a,e17.11)') & ' rho_sort(',nsort,') = ',rho_sort(nsort), & ' vol_sort(',nsort,') = ',vol_sort(nsort)*rho0r enddo endif ! layer properties, where layers defined by fixed number of points per layer nsort=0 do k=1,nk do n=1,nlayer(k) nsort=nsort+1 dzt_lay(k,2) = dzt_lay(k,2) + vol_sort(nsort) dzt_rho_lay(k,2) = dzt_rho_lay(k,2) + vol_sort(nsort)*rho_sort(nsort) enddo dzt_lay(k,2) = dzt_lay(k,2)*normalize(k) dzt_rho_lay(k,2) = dzt_rho_lay(k,2)*normalize(k) enddo deallocate(rho_sort) deallocate(vol_sort) do k=1,nk volume_lay(k,2) = dzt_lay(k,2)*wetarea(k) mass_lay(k,2) = dzt_rho_lay(k,2)*wetarea(k) if(volume_lay(k,2) > 0) then rho_lay(k,2) = mass_lay(k,2)/volume_lay(k,2) endif enddo do k=1,nk-1 dzw_lay(k,2) = 0.5*(dzt_lay(k,2)+dzt_lay(k+1,2)) enddo dzw_lay(nk,2) = 0.5*dzt_lay(nk,2) do k=1,nk thick_column(2) = thick_column(2) + dzt_lay(k,2) mass_column(2) = mass_column(2) + mass_lay(k,2) enddo ! compute derivatives across averaged layers. ! derivatives are defined at bottom of tracer cell, ! with deriv_lay(k=1) at bottom of the top-most cell, ! and deriv_lay(k=nk) at bottom of the bottom-most cell (and so set to zero). deriv_lay(:) = 0.0 do k=1,nk-1 if(dzw_lay(k,1) /= 0.0) then deriv_lay(k) = -(rho_lay(k+1,1)-rho_lay(k,1))/dzw_lay(k,1) endif enddo ! Compute diapycnal flux of sorted density. ! Flux is defined on the bottom face of the t-cell. ! Start at column top (k=1), specifying no flux condition. ! Bottom (k=nk) should be diagnosed to have zero flux for case ! where there is no buoyancy forcing. Finding a zero flux ! at bottom is a useful check on algorithm integrity. flux_lay(:) = 0.0 k=1 flux_lay(k) = (dzt_rho_lay(k,2)-dzt_rho_lay(k,1))*dtimer do k=2,nk flux_lay(k) = flux_lay(k-1) + (dzt_rho_lay(k,2)-dzt_rho_lay(k,1))*dtimer enddo ! compute effective diffusivity kappa_lay(:) = 0.0 numzero = 0 do k=1,nk if( abs(deriv_lay(k)) >= rho_grad_min .and. abs(deriv_lay(k)) <= rho_grad_max ) then kappa_lay(k) = -flux_lay(k)/deriv_lay(k) else numzero = numzero + 1 endif enddo if(smooth_kappa_sort>0) then do m=1,smooth_kappa_sort kappa_prev = 0.25*kappa_lay(1) do k=2,nk-2 tmp = kappa_lay(k) kappa_lay(k) = kappa_prev + 0.5*kappa_lay(k) + 0.25*kappa_lay(k+1) kappa_prev = 0.25*tmp enddo enddo endif if(debug_diagnose_mixingA) then write(stdoutunit,'(/a/)')'==========Debug diagnostics for kappa_sort calculation=========' write(stdoutunit,'(/a)')' Number of parcels within each density layer' do k=1,nk write(stdoutunit,'(a,i3,a,i7)') ' nlayer(', k,') = ',nlayer(k) enddo write(stdoutunit,'(/a)')' Thickness(m) of density layers' do k=1,nk difference = dzt_lay(k,2)-dzt_lay(k,1) write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11,a,i3,a,e17.11)') & ' dzt_lay(',k,',1)= ',dzt_lay(k,1), & ' dzt_lay(',k,',2)= ',dzt_lay(k,2), & ' diff(',k,')= ',difference enddo write(stdoutunit,'(/a)') ' Total thickness(m) of sorted column' write(stdoutunit,'(a,e17.11,a,e17.11,a,e17.11)') & ' thick_column(tau)= ' ,thick_column(tau) , & ' thick_column(taup1)= ',thick_column(taup1), & ' diff= ', thick_column(taup1)-thick_column(tau) write(stdoutunit,'(/a)')' Distance(m) between density cell centers ' do k=1,nk difference = dzw_lay(k,2)-dzw_lay(k,1) write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11,a,i3,a,e17.11)') & ' dzw_lay(',k,',1)= ',dzw_lay(k,1), & ' dzw_lay(',k,',2)= ',dzw_lay(k,2), & ' diff(',k,')= ',difference enddo write(stdoutunit,'(/a)') ' Mass (kg) of layers' do k=1,nk difference = mass_lay(k,2)-mass_lay(k,1) write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11,a,i3,a,e17.11)') & ' mass_lay(',k,',1)= ',mass_lay(k,1), & ' mass_lay(',k,',2)= ',mass_lay(k,2), & ' diff(',k,')= ',difference enddo write(stdoutunit,'(/a)') ' Total mass(kg) of sorted column' write(stdoutunit,'(a,e17.11,a,e17.11,a,e17.11)') & ' mass_column(tau)= ' ,mass_column(tau) , & ' mass_column(taup1)= ',mass_column(taup1), & ' diff= ', mass_column(taup1)-mass_column(tau) write(stdoutunit,'(/a)') ' Volume (m^3) of layers' do k=1,nk difference = volume_lay(k,2)-volume_lay(k,1) write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11,a,i3,a,e17.11)') & ' volume_lay(',k,',1)= ',volume_lay(k,1), & ' volume_lay(',k,',2)= ',volume_lay(k,2), & ' diff(',k,')= ',difference enddo write(stdoutunit,'(/a)') ' Density (kg/m^3) of layers' do k=1,nk difference = rho_lay(k,2)-rho_lay(k,1) write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11,a,i3,a,e17.11)') & ' rho_lay(',k,',1)= ',rho_lay(k,1), & ' rho_lay(',k,',2)= ',rho_lay(k,2), & ' diff(',k,')= ',difference enddo write(stdoutunit,'(/a)') ' Thickness weighted density dzt*rho (kg/m^2)' do k=1,nk difference = dzt_rho_lay(k,2)-dzt_rho_lay(k,1) write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11,a,i3,a,e17.11)') & ' dzt_rho_lay(',k,',1)= ',dzt_rho_lay(k,1), & ' dzt_rho_lay(',k,',2)= ',dzt_rho_lay(k,2), & ' diff(',k,')= ',difference enddo write(stdoutunit,'(/a)')' Terms contributing to effective diffusivity' do k=1,nk write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11)') & ' flux_lay(',k,') = ',flux_lay(k), & ' deriv_lay(',k,') = ',deriv_lay(k) enddo write (stdoutunit,'(/a)') ' Effective diapycnal diffusivity (m^2/sec)' do k=1,nk write(stdoutunit,'(a,i3,a,e17.11)')' kappa_sort(',k,') = ',kappa_lay(k) enddo write(stdoutunit,'(a)') '=============================================================' endif if (id_kappa_sort > 0) then used = send_data (id_kappa_sort, kappa_lay(:), Time%model_time) endif if (id_temp_sort > 0) then used = send_data (id_temp_sort, -rho_lay(:,1)/alpha, Time%model_time) endif end subroutine diagnose_kappa_sort ! NAME="diagnose_kappa_sort" !####################################################################### ! ! ! ! Routine to diagnose the amount of mixing between classes of a ! particular tracer. Temperature is used as default. ! ! Compute horizontal average of temp to define a stable profile. ! Evolution of this profile defines an effective diffusity. ! ! This diffusivity is different than the one diagnosed ! from the adiabatic sorting approach. The sorting approach is ! more relevant. The two approaches agree when there ! is zero baroclinicity, and the present simple scheme is ! useful ONLY to help debug the sorting routine. ! ! ! subroutine diagnose_kappa_simple(Time, Theta) type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), intent(in) :: Theta integer :: tau, taum1, taup1 real :: kappa_simple(nk) ! diagnosed dianeutral diffusivity (m^2/sec) real :: deriv_simple(nk) ! vertical density gradient centered at top of sorted tracer cell real :: flux_simple(nk) ! dianuetral tracer flux centered at top of sorted tracer cell real :: rho_simple(nk,3) ! density real, dimension(:,:), allocatable :: global_dat ! for global area real, dimension(:,:,:), allocatable :: global_tmask ! for global mask real, dimension(:,:,:,:), allocatable :: global_tracer ! for global tracer integer :: i, j, k real :: cellarea_r integer :: stdoutunit stdoutunit=stdout() if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (diagnose_kappa_simple): module needs initialization ') endif tau = Time%tau taum1 = Time%taum1 taup1 = Time%taup1 ! NOTE: minimum vertical density gradient rho_grad_min is necessary to ! avoid errors with truncation in the division by drho/dz when compute kappa. ! rho_grad_min corresponds roughly to the precision of the computation. ! Physically, with ! ! N^2 = -(g/rho0)(drho/dz) ! ! then rho_grad_min sets a minimum N^2 resolved. ! This corresponds to a frequency f=N/2pi. The typical ! period of inertial oscillations in the deep ocean is 6hrs ! (Pickard and Emery, page 56). In the upper ocean, it is ! 10-30 minutes. So to cover the majority of the ocean's stratification, ! we will want to set rho_grad_min to something smaller than 9e-6. ! ! To bin the effective diffusivity, it is also useful to have a max ! vertical density gradient. allocate (global_dat(ni,nj)) ; global_dat=0.0 call mpp_global_field(Dom%domain2d, Grd%dat, global_dat) allocate (global_tmask(ni,nj,nk)) ; global_tmask=0.0 call mpp_global_field(Dom%domain2d, Grd%tmask, global_tmask) ! horizontal tracer area if(debug_diagnose_mixingA) then write(stdoutunit,'(a)') ' ' write(stdoutunit,'(a,e17.11)')' in simple, total cross area of domain = ',Grd%tcellsurf write(stdoutunit,'(a,e17.11)')' in simple, total volume of domain = ',Grd%tcellsurf*Grd%zw(nk) endif cellarea_r = 1.0/(epsln + Grd%tcellsurf) ! horizontally average to get vertical density profile allocate (global_tracer(ni,nj,nk,3)) ; global_tracer=0.0 call mpp_global_field(Dom%domain2d, Theta%field(:,:,:,:), global_tracer) rho_simple(:,:) = 0.0 do k=1,nk do j=1,nj do i=1,ni if(global_tmask(i,j,k) == 1.0) then rho_simple(k,:) = -alpha*global_tracer(i,j,k,:)*global_dat(i,j) + rho_simple(k,:) endif enddo enddo rho_simple(k,:) = rho_simple(k,:)*cellarea_r enddo deallocate(global_tracer) deallocate(global_tmask) deallocate(global_dat) ! compute derivatives across averaged layers. ! derivatives are defined at bottom of tracer cell, ! with deriv_simple(k=1) at bottom of top-most cell, ! and deriv_simple(k=nk)=0 at bottom of bottom-most cell. deriv_simple(:) = 0.0 do k=1,nk-1 deriv_simple(k) = -(rho_simple(k+1,taup1)-rho_simple(k,taup1))*Grd%dzwr(k-1) enddo if(debug_diagnose_mixingA) then write(stdoutunit,'(/a)') 'density at three time levels' do k=1,nk write(stdoutunit,'(a,i3,a,e17.11,a,i3,a,e17.11,a,i3,a,e17.11)') & ' rho_simple(taum1,',k,')= ',rho_simple(k,taum1)+rho0, & ' rho_simple(tau(',k,')= ',rho_simple(k,tau)+rho0,& ' rho_simple(taup1(',k,')= ',rho_simple(k,taup1)+rho0 enddo write(stdoutunit,'(/a)') 'vertical derivative of density' do k=1,nk write(stdoutunit,'(a,i3,a,e17.11)')' deriv_simple(',k,') = ',deriv_simple(k) enddo endif ! compute effective diffusivity by diagnosing the ! diapycnal flux. The flux is defined on the bottom ! face of the t-cell. start integration from the ! ocean top and work down. flux_simple(:) = 0.0 k=1 flux_simple(k) = (rho_simple(k,taup1)-rho_simple(k,taum1))*Grd%dzt(k)*dtimer do k=2,nk flux_simple(k) = flux_simple(k-1) + (rho_simple(k,taup1)-rho_simple(k,taum1))*Grd%dzt(k)*dtimer enddo if(debug_diagnose_mixingA) then write(stdoutunit,'(/a)')'vertical flux across a layer' do k=1,nk write(stdoutunit,'(a,i3,a,e17.11)')' flux_simple(',k,') = ',flux_simple(k) enddo endif do k=1,nk if(abs(deriv_simple(k)) >= rho_grad_min .and. abs(deriv_simple(k)) <= rho_grad_max) then kappa_simple(k) = -flux_simple(k)/deriv_simple(k) else kappa_simple(k) = 0.0 endif enddo if (id_kappa_simple > 0) then used = send_data (id_kappa_simple, kappa_simple(:), Time%model_time) endif if(debug_diagnose_mixingA) then write (stdoutunit,'(/a)') ' Effective diffusivity (m^2/sec)' do k=1,nk write(stdoutunit,'(a,i3,a,e17.11)')' kappa_simple(',k,') = ',kappa_simple(k) enddo endif end subroutine diagnose_kappa_simple ! NAME="diagnose_kappa_simple" !####################################################################### ! ! ! ! Diagnose depth (m) of a potential density surface surface relative to ! the ocean surface at z=eta (not relative to z=0). ! ! Method uses linear interpolation to find the depth of a ! potential rho surface. ! ! Scheme currently does not forward (backwards) interpolate if ! rho surface lies within lowest (uppermost) grid cell. ! ! Diagnostic only makes sense when rho is monotonically ! increasing as go deeper in water column. ! ! Author: Harper.Simmons ! Zhi.Liang ! ! subroutine diagnose_depth_of_potrho(Time, Dens, Thickness) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens type(ocean_thickness_type), intent(in) :: Thickness real :: w1,w2 integer :: potrho_nk integer :: i, j, k, n real, dimension(isd:ied,jsd:jed,size(Dens%potrho_ref(:))) :: depth_of_potrho potrho_nk = size(Dens%potrho_ref(:)) if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (diagnose_depth_of_potrho): module needs initialization ') endif depth_of_potrho(:,:,:)=-10.0 do n=1,potrho_nk do j=jsc,jec do i=isc,iec kloop: do k=nk-1,1,-1 if( Dens%potrho(i,j,k) < Dens%potrho_ref(n)) then if(Dens%potrho_ref(n) < Dens%potrho(i,j,k+1)) then if(Grd%tmask(i,j,k+1) > 0) then W1= Dens%potrho_ref(n) - Dens%potrho(i,j,k) W2= Dens%potrho(i,j,k+1) - Dens%potrho_ref(n) depth_of_potrho(i,j,n) = ( Thickness%depth_zt(i,j,k+1)*W1 & +Thickness%depth_zt(i,j,k) *W2) & /(W1 + W2 + epsln) exit kloop endif endif endif enddo kloop enddo enddo enddo used = send_data (id_depth_of_potrho, depth_of_potrho(:,:,:), & Time%model_time, & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=potrho_nk) end subroutine diagnose_depth_of_potrho ! NAME="diagnose_depth_of_potrho" !####################################################################### ! ! ! ! Diagnose depth (m) of a potential temperature surface relative to ! the ocean surface at z=eta (not relative to z=0). ! ! Method uses linear interpolation to find the depth of a ! potential temp surface. ! ! Scheme currently does not forward (backwards) interpolate if ! theta surface lies within lowest (uppermost) grid cell. ! ! Diagnostic only makes sense when theta is monotonically ! decreasing as go deeper in water column. ! ! Based on "diagnose_depth_of_potrho" by Harper.Simmons ! ! Author: Stephen.Griffies ! Zhi.Liang ! ! subroutine diagnose_depth_of_theta(Time, Dens, Thickness, T_prog) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: T_prog(:) real :: w1,w2 integer :: theta_nk integer :: i, j, k, n, tau real, dimension(isd:ied,jsd:jed,size(Dens%theta_ref(:))) :: depth_of_theta theta_nk = size(Dens%theta_ref(:)) tau = Time%tau if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (diagnose_depth_of_theta): module needs initialization ') endif depth_of_theta(:,:,:)=-10.0 do n=1,theta_nk do j=jsc,jec do i=isc,iec kloop: do k=nk-1,1,-1 if( Dens%theta_ref(n) < T_prog(index_temp)%field(i,j,k,tau) ) then if(T_prog(index_temp)%field(i,j,k+1,tau) < Dens%theta_ref(n)) then if(Grd%tmask(i,j,k+1) > 0) then W1= T_prog(index_temp)%field(i,j,k,tau) - Dens%theta_ref(n) W2= Dens%theta_ref(n) - T_prog(index_temp)%field(i,j,k+1,tau) depth_of_theta(i,j,n) = ( Thickness%depth_zt(i,j,k+1)*W1 & +Thickness%depth_zt(i,j,k) *W2 ) & /(W1 + W2 + epsln) exit kloop endif endif endif enddo kloop enddo enddo enddo used = send_data (id_depth_of_theta, depth_of_theta(:,:,:), & Time%model_time, & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=theta_nk) end subroutine diagnose_depth_of_theta ! NAME="diagnose_depth_of_theta" !####################################################################### ! ! ! Diagnose tracer concentration on potential density surface. ! Method based on diagnose_depth_of_potrho diagnostic. ! ! Author: Stephen.Griffies ! ! Updated Oct 2009 to be more vectorized ! ! ! subroutine diagnose_tracer_on_rho(Time, Dens, Tracer, ntracer) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens type(ocean_prog_tracer_type), intent(in) :: Tracer integer, intent(in) :: ntracer real :: W1, W2 integer :: potrho_nk integer :: i, j, k, k_rho, tau real, dimension(isd:ied,jsd:jed,size(Dens%potrho_ref(:))) :: tracer_on_rho real, dimension(jsc:jec) :: rho_minj,rho_maxj real :: rho_min,rho_max potrho_nk = size(Dens%potrho_ref(:)) tau = Time%tau if (.not. module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (diagnose_tracer_on_rho): module needs initialization') endif tracer_on_rho(:,:,:)=0.0 ! for (i,j) points with potrho_ref < potrho(k=1), keep tracer_on_rho=0 ! for (i,j) points with potrho_ref > potrho(k=kmt), keep tracer_on_rho=0 ! these assumptions mean there is no need to specially handle the endpoints, ! since the initial value for tracer_on_rho is 0. ! interpolate tracer field onto rho-surface do k = 1,nk-1 do j = jsc,jec rho_maxj(j) = maxval(Dens%potrho(isc:iec,j,k+1),mask=Grd%tmask(isc:iec,j,k+1)==1.) rho_minj(j) = minval(Dens%potrho(isc:iec,j,k),mask=Grd%tmask(isc:iec,j,k)==1.) enddo rho_max = maxval(rho_maxj) rho_min = minval(rho_minj) if (rho_max == -huge(rho_max)) exit ! only rock below this level do k_rho = 1,potrho_nk if (rho_max < Dens%potrho_ref(k_rho)) cycle if (rho_min > Dens%potrho_ref(k_rho)) cycle do j = jsc,jec if (rho_maxj(j) < Dens%potrho_ref(k_rho)) cycle if (rho_minj(j) > Dens%potrho_ref(k_rho)) cycle do i=isc,iec if( Dens%potrho_ref(k_rho) > Dens%potrho(i,j,k) ) then if( Dens%potrho_ref(k_rho) <= Dens%potrho(i,j,k+1)) then W1= Dens%potrho_ref(k_rho)- Dens%potrho(i,j,k) W2= Dens%potrho(i,j,k+1) - Dens%potrho_ref(k_rho) tracer_on_rho(i,j,k_rho) = ( Tracer%field(i,j,k+1,tau)*W1 & +Tracer%field(i,j,k,tau) *W2) & /(W1 + W2 + epsln) endif endif enddo enddo enddo enddo ! ensure masking is applied to interpolated field do k_rho = 1,potrho_nk do j = jsc,jec do i = isc,iec tracer_on_rho(i,j,k_rho) = tracer_on_rho(i,j,k_rho)*Grd%tmask(i,j,1) enddo enddo enddo used = send_data (id_tracer_on_rho(ntracer), tracer_on_rho(:,:,:), & Time%model_time, & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=potrho_nk) end subroutine diagnose_tracer_on_rho ! NAME="diagnose_tracer_on_rho" !####################################################################### ! ! ! Diagnose tracer concentration * dz/drho on potential density surface. ! This product, when integrated over dx*dy*drho, will yield the same ! total tracer (to within roundoff) as the usual tracer concentration ! integrated over dx*dy*dz. ! ! compute abs(dz/drho)==dz/drho in order to have tracer_zrho_on_rho ! with same sign as tracer. ! ! Method based on diagnose_tracer_on_rho diagnostic. ! ! Author: Stephen.Griffies ! Updated Oct 2009 to be more vectorized ! ! ! subroutine diagnose_tracer_zrho_on_rho(Time, Dens, Thickness, Tracer, ntracer) type(ocean_time_type), intent(in) :: Time type(ocean_density_type), intent(in) :: Dens type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Tracer integer, intent(in) :: ntracer real :: W1,W2 real :: zrho integer :: potrho_nk integer :: i, j, k, n, tau real, dimension(isd:ied,jsd:jed,size(Dens%potrho_ref(:))) :: tracer_zrho_on_rho potrho_nk = size(Dens%potrho_ref(:)) tau = Time%tau if (.not. module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (diagnose_tracer_zrho_on_rho): module needs initialization') endif tracer_zrho_on_rho(:,:,:)=0.0 ! for (i,j) points with potrho_ref < potrho(k=1), keep tracer_zrho_on_rho=0 ! for (i,j) points with potrho_ref > potrho(k=kmt), keep tracer_zrho_on_rho=0 ! these assumptions mean there is no need to specially handle the endpoints, ! since the initial value for tracer_on_rho is 0. ! interpolate field onto rho-surface do n=1,potrho_nk do k=nk-1,1,-1 do j=jsc,jec do i=isc,iec if( Dens%potrho(i,j,k) < Dens%potrho_ref(n)) then if(Dens%potrho_ref(n) < Dens%potrho(i,j,k+1)) then if(Grd%tmask(i,j,k+1) > 0) then zrho = Thickness%dzwt(i,j,k)/(-Dens%potrho(i,j,k)+Dens%potrho(i,j,k+1)+epsln) W1= Dens%potrho_ref(n) - Dens%potrho(i,j,k) W2= Dens%potrho(i,j,k+1) - Dens%potrho_ref(n) tracer_zrho_on_rho(i,j,n) = zrho* & ( Tracer%field(i,j,k+1,tau)*W1 & +Tracer%field(i,j,k,tau) *W2) & /(W1 + W2 + epsln) endif endif endif enddo enddo enddo enddo used = send_data (id_tracer_zrho_on_rho(ntracer), tracer_zrho_on_rho(:,:,:), & Time%model_time, & is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=potrho_nk) end subroutine diagnose_tracer_zrho_on_rho ! NAME="diagnose_tracer_zrho_on_rho" !####################################################################### ! ! ! ! Calculate the mixed layer depth and potential density at mixed layer base ! according to depth at which buoyancy is greater than buoyancy_crit ! relative to the surface. Compute the buoyancy using potential ! density, rather than the insitu density, since we aim for this ! diagnostic to be comparable to diagnostics from isopcynal models. ! ! Note that the mixed layer depth is taken with respect to the ocean surface, ! and so the mixed layer depth is always positive. That is, the mld is here ! defined as a thickness of water. ! ! ! subroutine calc_potrho_mixed_layer (Thickness, Dens, potrho_mix_depth, potrho_mix_base) type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens real, dimension(isd:,jsd:), intent(out), optional :: potrho_mix_depth real, dimension(isd:,jsd:), intent(out), optional :: potrho_mix_base integer :: i, j, k real :: depth, crit real, dimension(isd:ied,jsd:jed) :: potrho_mix_base_use if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (calc_potrho_mixed_layer): module needs initialization ') endif crit = buoyancy_crit*rho0/grav potrho_mix_base_use = 0.0 do j=jsc,jec do i=isc,iec potrho_mix_base_use(i,j) = Dens%potrho(i,j,1) + crit enddo enddo if (present(potrho_mix_base)) then potrho_mix_base = potrho_mix_base_use endif if (present(potrho_mix_depth)) then potrho_mix_depth = 0.0 do j=jsc,jec do i=isc,iec depth=max(0.0,Thickness%dzwt(i,j,0)) do k=2,Grd%kmt(i,j) if( Dens%potrho(i,j,k) < potrho_mix_base_use(i,j)) then depth=depth+Thickness%dzwt(i,j,k) else potrho_mix_depth(i,j) = depth + Thickness%dzwt(i,j,k) & *(potrho_mix_base_use(i,j)-Dens%potrho(i,j,k-1))& /(Dens%potrho(i,j,k) -Dens%potrho(i,j,k-1)) exit endif enddo enddo enddo endif end subroutine calc_potrho_mixed_layer ! NAME="calc_potrho_mixed_layer" !####################################################################### ! ! ! ! Determine mixed layer depth and potential density at mixed layer base ! according to depth at which buoyancy is greater than buoyancy_crit ! relative to the surface. ! Call calc_potrho_mixed_layer to calculate the quantities. ! ! ! subroutine potrho_mixed_layer (Time, Thickness, Dens) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens type(time_type) :: next_time real :: potrho_mix_depth(isd:ied,jsd:jed) real :: potrho_mix_base(isd:ied,jsd:jed) if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (potrho_mixed_layer): module needs initialization ') endif next_time = increment_time(Time%model_time, int(dtts), 0) if (need_data(id_potrho_mix_depth,next_time) .and. id_potrho_mix_base > 0) then call calc_potrho_mixed_layer (Thickness, Dens, & potrho_mix_depth = potrho_mix_depth, potrho_mix_base = potrho_mix_base) elseif (need_data(id_potrho_mix_depth,next_time)) then call calc_potrho_mixed_layer (Thickness, Dens, & potrho_mix_depth = potrho_mix_depth) elseif (id_potrho_mix_base > 0) then call calc_potrho_mixed_layer (Thickness, Dens, & potrho_mix_base = potrho_mix_base) endif call diagnose_2d(Time, Grd, id_potrho_mix_depth, potrho_mix_depth(:,:)) call diagnose_2d(Time, Grd, id_potrho_mix_base, potrho_mix_base(:,:)) end subroutine potrho_mixed_layer ! NAME="potrho_mixed_layer" !####################################################################### ! ! ! ! Send total liquid seawater mass to diagnostic manager. ! ! subroutine send_total_mass(Time, Thickness) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness integer :: tau real :: tmass if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (send_total_mass): module needs initialization ') endif tau = Time%tau tmass = total_mass(Thickness,tau) tmass = tmass used = send_data (id_mass_seawater, tmass, Time%model_time) end subroutine send_total_mass ! NAME="send_total_mass" !####################################################################### ! ! ! ! Send total liquid seawater mass to diagnostic manager. ! ! subroutine send_total_volume(Time, Thickness) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness integer :: tau real :: tvolume tau = Time%tau if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (send_total_volume): module needs initialization ') endif tvolume = total_volume(Thickness, tau) tvolume = tvolume used = send_data (id_volume_seawater, tvolume, Time%model_time) end subroutine send_total_volume ! NAME="send_total_volume" !####################################################################### ! ! ! ! Send total tracer to diagnostic manager. ! ! subroutine send_total_tracer(Time, Thickness, Tracer, ntracer) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Tracer integer, intent(in) :: ntracer integer :: tau real :: ttracer if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (send_total_tracer): module needs initialization ') endif tau = Time%tau ttracer = total_tracer(Tracer,Thickness,tau) if(Tracer%name=='temp') then ttracer = ttracer*1e-25 elseif(Tracer%name(1:3)=='age') then ttracer = ttracer else ttracer = ttracer*1e-18 endif used = send_data (id_prog_total(ntracer), ttracer, Time%model_time) end subroutine send_total_tracer ! NAME="send_total_tracer" !####################################################################### ! ! ! ! Send global averaged tracer to diagnostic manager. ! ! subroutine send_global_ave_tracer(Time, Thickness, Tracer, ntracer) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Tracer integer, intent(in) :: ntracer real :: ave_tracer, totaltracer, totalmass integer :: tau if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (send_global_ave_tracer): module needs initialization ') endif tau=Time%tau totaltracer = total_tracer(Tracer,Thickness,tau) totalmass = total_mass(Thickness,tau) ave_tracer = totaltracer/totalmass/Tracer%conversion used = send_data (id_prog_global_ave(ntracer), ave_tracer, Time%model_time) end subroutine send_global_ave_tracer ! NAME="send_global_ave_tracer" !####################################################################### ! ! ! ! Send global averaged pressure to diagnostic manager. ! ! subroutine send_global_ave_pressure(Time, Thickness, Dens) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens real :: total_press, press_ave, totalmass real, dimension(isd:ied,jsd:jed) :: press_k integer :: tau, k if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (send_global_ave_press): module needs initialization') endif tau=Time%tau total_press = 0.0 press_k(:,:) = 0.0 do k=1,nk press_k(:,:) = press_k(:,:) + Grd%tmask(:,:,k)*Grd%dat(:,:) & *Thickness%rho_dzt(:,:,k,tau)*Dens%pressure_at_depth(:,:,k) enddo if(have_obc) press_k(:,:) = press_k(:,:)*Grd%obc_tmask(:,:) total_press = mpp_global_sum(Dom%domain2d, press_k(:,:), global_sum_flag) totalmass = total_mass(Thickness,tau) press_ave = total_press/totalmass used = send_data (id_press_ave, press_ave, Time%model_time) end subroutine send_global_ave_pressure ! NAME="send_global_ave_pressure" !####################################################################### ! ! ! ! Send global averaged surface tracer to diagnostic manager. ! Note the presence of a rho_dzt weighting here... ! ! subroutine send_surface_ave_tracer(Time, Thickness, Tracer, ntracer) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_prog_tracer_type), intent(in) :: Tracer integer, intent(in) :: ntracer real :: ave_tracer, totaltracer, totalmass integer :: tau if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (send_surface_ave_tracer): module needs initialization ') endif tau=Time%tau totaltracer = klevel_total_tracer(Tracer,Thickness,tau,1) totalmass = klevel_total_mass(Thickness,tau,1) ave_tracer = totaltracer/totalmass/Tracer%conversion used = send_data (id_prog_surface_ave(ntracer), ave_tracer, Time%model_time) end subroutine send_surface_ave_tracer ! NAME="send_surface_ave_tracer" !####################################################################### ! ! ! ! Send global area averaged surface tracer to diagnostic manager. ! Note the weigthing is just area, with no thickness nor density. ! ! subroutine send_surface_area_ave_tracer(Time, Tracer, ntracer) type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), intent(in) :: Tracer integer, intent(in) :: ntracer real, dimension(isd:ied,jsd:jed) :: tmp real :: ave_tracer integer :: i,j,tau if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error from ocean_tracer_diag (send_surface_area_ave_tracer): module needs initialization') endif tau=Time%tau tmp(:,:) = 0.0 do j=jsc,jec do i=isc,iec tmp(i,j) = tmp(i,j) + Tracer%field(i,j,1,tau)*area_t(i,j) enddo enddo ave_tracer = mpp_global_sum(Dom%domain2d, tmp(:,:), global_sum_flag) ave_tracer = ave_tracer*cellarea_r used = send_data (id_prog_surface_area_ave(ntracer), ave_tracer, Time%model_time) end subroutine send_surface_area_ave_tracer ! NAME="send_surface_area_ave_tracer" !####################################################################### ! ! ! ! ! Compute the global and k-level tracer variance. ! ! ! subroutine send_tracer_variance(Time, Tracer, Thickness, n) type(ocean_time_type), intent(in) :: Time type(ocean_prog_tracer_type), intent(in) :: Tracer type(ocean_thickness_type), intent(in) :: Thickness integer, intent(in) :: n integer :: i,j,k,tau real :: global_variance, tmass, ttracer, ttracer2 real, dimension(nk) :: k_variance, k_tmass, k_ttracer, k_ttracer2 tau = Time%tau if (id_global_variance(n) > 0 .or. id_k_variance(n) > 0) then if (use_blobs) then wrk1(:,:,:) = 0.0 wrk2(:,:,:) = 0.0 wrk3(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec wrk1(i,j,k) = Grd%dat(i,j)*Thickness%rho_dztT(i,j,k,tau)*Grd%tmask(i,j,k) !Global Mass wrk2(i,j,k) = wrk1(i,j,k)*Tracer%fieldT(i,j,k)*Tracer%fieldT(i,j,k) wrk3(i,j,k) = wrk1(i,j,k)*Tracer%fieldT(i,j,k) enddo enddo enddo else!not use_blobs wrk1(:,:,:) = 0.0 wrk2(:,:,:) = 0.0 wrk3(:,:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec wrk1(i,j,k) = Grd%dat(i,j)*Thickness%rho_dzt(i,j,k,tau)*Grd%tmask(i,j,k) wrk2(i,j,k) = wrk1(i,j,k)*Tracer%field(i,j,k,tau)*Tracer%field(i,j,k,tau) wrk3(i,j,k) = wrk1(i,j,k)*Tracer%field(i,j,k,tau) enddo enddo enddo endif ! endif for blobs if (id_global_variance(n) > 0) then tmass = mpp_global_sum(Dom%domain2d, wrk1(:,:,:)) ttracer2 = mpp_global_sum(Dom%domain2d, wrk2(:,:,:)) ttracer = mpp_global_sum(Dom%domain2d, wrk3(:,:,:)) if(tmass > 0.0) then global_variance = ttracer2/tmass - ( ttracer/tmass )**2 else global_variance = 0.0 endif used = send_data (id_global_variance(n), global_variance, Time%model_time) endif if (id_k_variance(n) > 0) then do k=1,nk k_tmass(k) = mpp_global_sum(Dom%domain2d, wrk1(:,:,k)) k_ttracer2(k) = mpp_global_sum(Dom%domain2d, wrk2(:,:,k)) k_ttracer(k) = mpp_global_sum(Dom%domain2d, wrk3(:,:,k)) if(k_tmass(k) > 0.0) then k_variance(k) = k_ttracer2(k)/k_tmass(k) - ( k_ttracer(k)/k_tmass(k) )**2 else k_variance(k) = 0.0 endif enddo used = send_data( id_k_variance(n), k_variance(:), Time%model_time ) endif endif end subroutine send_tracer_variance ! NAME="send_tracer_variance" !####################################################################### ! ! ! ! Diagnose contribution to global mean sea level evolution arising ! from a 3d MOM flux computed from a parameterization. ! ! fluxes are assumed to have the following dimensions: ! flux_x = (dy*dz)*diffusivity*rho*tracer_derivative_x ! flux_y = (dx*dz)*diffusivity*rho*tracer_derivative_y ! flux_z = diffusivity*rho*tracer_derivative_z ! ! Subroutine history: ! Jan2012 version 1.0: Stephen.Griffies ! ! ! subroutine diagnose_eta_tend_3dflux (Time, Thickness, Dens,& flux_x_temp, flux_y_temp, flux_z_temp, & flux_x_salt, flux_y_salt, flux_z_salt, & eta_tend, eta_tend_glob) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness type(ocean_density_type), intent(in) :: Dens real, dimension(isd:,jsd:,:), intent(in) :: flux_x_temp real, dimension(isd:,jsd:,:), intent(in) :: flux_y_temp real, dimension(isd:,jsd:,:), intent(in) :: flux_z_temp real, dimension(isd:,jsd:,:), intent(in) :: flux_x_salt real, dimension(isd:,jsd:,:), intent(in) :: flux_y_salt real, dimension(isd:,jsd:,:), intent(in) :: flux_z_salt real, dimension(isd:,jsd:), intent(inout) :: eta_tend real, intent(inout) :: eta_tend_glob integer :: i, j, k, kp1, tau if (.not.module_is_initialized) then call mpp_error(FATAL, & '==>Error ocean_tracer_util_mod (diagnose_eta_tend_3dflux): module needs initialization ') endif tau = Time%tau wrk1(:,:,:) = 0.0 ! nu_theta = alpha/rho = -drhodT/rho**2 wrk2(:,:,:) = 0.0 ! nu_salinity = -beta/rho = -drhodS/rho**2 do k=1,nk do j=jsd,jed do i=isd,ied wrk1(i,j,k) = -Grd%tmask(i,j,k)*Dens%drhodT(i,j,k)/(epsln+Dens%rho(i,j,k,tau)**2) wrk2(i,j,k) = -Grd%tmask(i,j,k)*Dens%drhodS(i,j,k)/(epsln+Dens%rho(i,j,k,tau)**2) enddo enddo enddo ! compute partial_x(nu_theta) and partial_x(nu_salinity) at east face of T-cell ! compute partial_y(nu_theta) and partial_y(nu_salinity) at north face of T-cell wrk1_v(:,:,:,1) = 0.0 ! partial_x(nu_theta) in 1/(deg C * rho * meter) wrk1_v(:,:,:,2) = 0.0 ! partial_x(nu_salinity) in 1/(ppt * rho * meter) wrk2_v(:,:,:,1) = 0.0 ! partial_y(nu_theta) in 1/(deg C * rho * meter) wrk2_v(:,:,:,2) = 0.0 ! partial_y(nu_salinity) in 1/(ppt * rho * meter) do k=1,nk wrk3(:,:,1) = wrk1(:,:,k) wrk4(:,:,1) = wrk2(:,:,k) wrk1_v(:,:,k,1) = FDX_T(wrk3(:,:,1))*FMX(Grd%tmask(:,:,k)) wrk1_v(:,:,k,2) = FDX_T(wrk4(:,:,1))*FMX(Grd%tmask(:,:,k)) wrk2_v(:,:,k,1) = FDY_T(wrk3(:,:,1))*FMY(Grd%tmask(:,:,k)) wrk2_v(:,:,k,2) = FDY_T(wrk4(:,:,1))*FMY(Grd%tmask(:,:,k)) enddo ! compute partial_z(nu_theta) and partial_z(nu_salinity) at bottom of T-cell wrk3_v(:,:,:,1) = 0.0 ! partial_z(nu_theta) in 1/(deg C * rho * meter) wrk3_v(:,:,:,2) = 0.0 ! partial_z(nu_salinity) in 1/(ppt * rho * meter) do k=1,nk-1 kp1=k+1 do j=jsc,jec do i=isc,iec wrk3_v(i,j,k,1) = Grd%tmask(i,j,kp1)*(wrk1(i,j,k)-wrk1(i,j,kp1))/Thickness%dzwt(i,j,k) wrk3_v(i,j,k,2) = Grd%tmask(i,j,kp1)*(wrk2(i,j,k)-wrk2(i,j,kp1))/Thickness%dzwt(i,j,k) enddo enddo enddo ! horizontal flux components have a dzt factor already ! included (thickness weighted fluxes). ! dytr and dxtr multipliers account for extra length ! factors appearing in the horizontal flux components. ! the vertical flux component is multiplied by dzt ! to anticipate the vertical integral, and to bring ! it into consistent dimensions with the horizontal flux ! components. wrk3(:,:,:) = 0.0 do k=1,nk kp1=min(nk,k+1) do j=jsc,jec do i=isc,iec wrk3(i,j,k) = Grd%tmask(i+1,j,k)*Grd%dytr(i,j) & *(flux_x_temp(i,j,k)*wrk1_v(i,j,k,1) + flux_x_salt(i,j,k)*wrk1_v(i,j,k,2)) & + Grd%tmask(i,j+1,k)*Grd%dxtr(i,j) & *(flux_y_temp(i,j,k)*wrk2_v(i,j,k,1) + flux_y_salt(i,j,k)*wrk2_v(i,j,k,2)) & + Grd%tmask(i,j,kp1)*Thickness%dzt(i,j,k) & *(flux_z_temp(i,j,k)*wrk3_v(i,j,k,1) + flux_z_salt(i,j,k)*wrk3_v(i,j,k,2)) enddo enddo enddo ! sea level tendency via a vertical integral. ! minus sign switches from MOM convention of an "upgradient" ! flux, to the desired downgradient form. eta_tend(:,:) = 0.0 do k=1,nk do j=jsc,jec do i=isc,iec eta_tend(i,j) = eta_tend(i,j) - wrk3(i,j,k) enddo enddo enddo ! global integral wrk1_2d(:,:) = Grd%tmask(:,:,1)*Grd%dat(:,:)*eta_tend(:,:) eta_tend_glob = mpp_global_sum(Dom%domain2d, wrk1_2d(:,:), global_sum_flag)*cellarea_r end subroutine diagnose_eta_tend_3dflux ! NAME="diagnose_eta_tend_3dflux" end module ocean_tracer_diag_mod