!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module ocean_tracer_diag_mod ! ! S.M. Griffies ! ! ! R.C. Pacanowski ! ! ! ! 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. ! ! ! Days over which time average the thickness weighted density before taking its ! time tendency for use in computing the effective diapycnal diffusivity. ! ! ! ! The realistic EOS used in mom4 requires salinity to ! use the Practical Salinity Scale (pss). This scale is ! also known as the Practical Salinity Unit (psu). ! Additionally, ocean measurements use the psu scale ! Hence, mom4 interprets its salinity as 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 mom4 carries as its ! prognostic salinity field. ! ! 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. ! ! ! ! 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: grav, 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: 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 use ocean_parameters_mod, only: TWO_LEVEL, THREE_LEVEL, missing_value use ocean_parameters_mod, only: rho0, rho0r, sec_in_yr, sec_in_yr_r 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_util_mod, only: write_timestamp use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk1_2d, wrk2_2d, wrk1_v, wrk2_v , wrk3_v 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 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_mixed_layer_depth 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 1.1.2.2.2.34.28.4.44.1.32.1.4.2 2009/12/01 21:51:38 smg Exp $' character (len=128) :: tagname = & '$Name: mom4p1_pubrel_dec2009_nnz $' ! for tracer conservation: need to know num_prog_tracers to allocate real, dimension(:,:,:), allocatable :: tracer_flux real, dimension(:,:,:), allocatable :: tracer_source real, dimension(:,:,:), allocatable :: tracer_eta_smooth real, dimension(:,:,:), allocatable :: tracer_pbot_smooth real, dimension(:,:) , allocatable :: material_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_global_ave integer :: id_mass_seawater=-1 integer :: id_volume_seawater=-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 ! 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 ! for potential density mixed layer integer :: id_potrho_mix_depth=-1 integer :: id_potrho_mix_base=-1 real :: buoyancy_crit=0.0003 ! (m^2/sec) critical buoyancy difference relative to surface ! for salt budgets real :: psu2ppt=1.004867 ! conversion from mom4's salinity in 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. public ocean_tracer_diag_init public ocean_tracer_diagnostics public calc_mixed_layer_depth public calc_potrho_mixed_layer private mixed_layer_depth 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_surface_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 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, T_prog, T_diag, Dens, & ver_coordinate_class, 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_prog_tracer_type), intent(in) :: T_prog(:) type(ocean_diag_tracer_type), intent(in) :: T_diag(:) type(ocean_density_type), intent(in) :: Dens integer, intent(in) :: ver_coordinate_class logical, intent(in) :: obc integer :: n, ioun, io_status, ierr 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) ioun = open_namelist_file() read(ioun, ocean_tracer_diag_nml, iostat=io_status) write (stdlogunit, ocean_tracer_diag_nml) write (stdoutunit,'(/)') write (stdoutunit, ocean_tracer_diag_nml) ierr = check_nml_error(io_status,'ocean_tracer_diag_nml') call close_file(ioun) 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 MOM4_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 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_global_ave(num_prog_tracers) ) id_prog_total= -1 id_prog_surface_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 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=(/0.0,1.e3/)) 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') 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/)) 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/)) 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/)) 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/)) 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 = 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_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/)) ! define clock integers 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) 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, 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(in) :: 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) :: 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)) 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%pressure_at_depth(isd:ied,jsd:jed,:), Time%model_time) endif call mpp_clock_end(id_mixed_layer_depth) call mpp_clock_begin(id_potrho_mixed_layer) call potrho_mixed_layer(Time, Thickness, Dens) call mpp_clock_end(id_potrho_mixed_layer) ! 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) ! 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 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/s2). ! ! 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, model_time, 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 type(time_type), intent(in) :: model_time real, dimension(isd:,jsd:), intent(out) :: 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. ! ! ! subroutine mixed_layer_depth(Thickness, salinity, theta, rho, pressure, model_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) :: pressure type(time_type), intent(in) :: model_time 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, model_time, hmxl) if(id_mld > 0) then used = send_data (id_mld, hmxl(:,:), model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif if(id_mld_sq > 0) then used = send_data (id_mld_sq, hmxl(:,:)**2, model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif end subroutine mixed_layer_depth ! NAME="mixed_layer_depth" !####################################################################### ! ! ! ! ! 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) :: 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 :: dmasstracer(0:nk) real :: eta_t_max, eta_t_min, eta_t_avg real :: anompbot_max, anompbot_min real :: frazil_heat real :: tracer_th_tend(0:nk), tracer_src(0:nk) real :: ttracer real :: mass_total, mass_source, mass_change real :: mass_taup1, mass_taup1_r, mass_taup1_keq1 real :: mass_eta_smooth, mass_pbot_smooth real :: conversion 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) 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 mass_total = total_mass(Thickness, tau) ! 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 do k=1,nk tmp(:,:) = Grd%tmask(:,:,k)*area_t(:,:) & *(Thickness%rho_dzt(:,:,k,taup1)-Thickness%rho_dzt(:,:,k,taum1)) mass_change = mass_change + mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) enddo ! mass input to ocean from sources do k=1,nk tmp(:,:) = dtime*area_t(:,:)*Thickness%mass_source(:,:,k) mass_source = mass_source + mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) enddo ! 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 ! (smg: may need to be modified for general vertical coordinates) 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' write (stdoutunit,'(a,es24.17,a)') ' Mass change = ',& mass_change,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of ocean tracer cells (tau) = ',& mass_total,' kg' 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 ice melt input = ',& melt_input,' kg' write (stdoutunit,'(a,es24.17,a)') ' Mass of river runoff water input = ',& runoff_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) 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(Time, T_prog(n), tmp, n, send_out =.true.) 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(:,:) = T_diag(index_frazil)%field(:,:,1)*area_t(:,:)/frazil_factor frazil_heat = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tracer_input_total = tracer_input_total + frazil_heat endif do k=1,nk area_k(:,:) = Grd%dat(:,:)*Grd%tmask(:,:,k) if(have_obc) area_k(:,:) = area_k(:,:)*Grd%obc_tmask(:,:) tmp(:,:) = area_k(:,:)*T_prog(n)%conversion* & (Thickness%rho_dzt(:,:,k,taup1)*T_prog(n)%field(:,:,k,taup1) & -Thickness%rho_dzt(:,:,k,taum1)*T_prog(n)%field(:,:,k,taum1)) dmasstracer(k) = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = area_k(:,:)*conversion*dtime*T_prog(n)%th_tendency(:,:,k) if(have_obc) tmp(:,:) = tmp(:,:)*Grd%obc_tmask(:,:) tracer_th_tend(k) = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = area_k(:,:)*conversion*dtime*T_prog(n)%source(:,:,k)*Thickness%rho_dzt(:,:,k,taup1) if(have_obc) tmp(:,:) = tmp(:,:)*Grd%obc_tmask(:,:) tracer_src(k) = mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) enddo ! 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(1) = tracer_th_tend(1) -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(1) = tracer_th_tend(1) - tracer_input_stf endif do k=1,nk ! sum integrals over all levels tracer_th_tend(0) = tracer_th_tend(0) + tracer_th_tend(k) tracer_src(0) = tracer_src(0) + tracer_src(k) dmasstracer(0) = dmasstracer(0) + dmasstracer(k) enddo ttracer = total_tracer(T_prog(n),Thickness, tau) tracer_error = dmasstracer(0)-tracer_input_total-tracer_th_tend(0)-tracer_src(0) 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' 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(0),' J' write (stdoutunit,'(a,es24.17,a)') ' Heat input via sources in th_tendency, or errors = ',& tracer_th_tend(0),' 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(0),' 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(0)-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' 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(0)*mass_taup1_r,' yr' write (stdoutunit,'(a,es24.17,a)') ' Age input via sources in th_tendency, or errors = ',& tracer_th_tend(0)*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(0)*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' 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(0),' kg' write (stdoutunit,'(a,es24.17,a)') ' Tracer input via sources in th_tendency, or errors = ',& tracer_th_tend(0),' 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(0),' 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(0),' 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 do k=1,nk tracer_k(:,:) = Tracer%conversion*Grd%tmask(:,:,k)*Grd%dat(:,:) & *Thickness%rho_dzt(:,:,k,index)*Tracer%field(:,:,k,index) if(have_obc) tracer_k(:,:) = tracer_k(:,:)*Grd%obc_tmask(:,:) total_tracer = total_tracer + mpp_global_sum(Dom%domain2d,tracer_k(:,:), global_sum_flag) enddo 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 tracer_k(:,:) = Tracer%conversion*Grd%tmask(:,:,kindex)*Grd%dat(:,:) & *Thickness%rho_dzt(:,:,kindex,tindex)*Tracer%field(:,:,kindex,tindex) 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 do k=1,nk mass_t(:,:) = Thickness%rho_dzt(:,:,k,index)*Grd%dat(:,:)*Grd%tmask(:,:,k) if(have_obc) mass_t(:,:) = mass_t(:,:)*Grd%obc_tmask(:,:) total_mass = total_mass + mpp_global_sum(Dom%domain2d,mass_t(:,:), global_sum_flag) enddo end function total_mass ! NAME="total_mass" !####################################################################### ! ! ! ! Compute total ocean tracer cell volume. ! ! function total_volume (Thickness) type(ocean_thickness_type), intent(in) :: Thickness 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 do k=1,nk volume_t(:,:) = Thickness%dzt(:,:,k)*Grd%dat(:,:)*Grd%tmask(:,:,k) if(have_obc) volume_t(:,:) = volume_t(:,:)*Grd%obc_tmask(:,:) total_volume = total_volume + mpp_global_sum(Dom%domain2d, volume_t(:,:), global_sum_flag) enddo 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 mass_t(:,:) = Thickness%rho_dzt(:,:,kindex,tindex)*Grd%dat(:,:)*Grd%tmask(:,:,kindex) 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) :: tmp real, dimension(isd:ied,jsd:jed) :: mass_t 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) 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(:,:) tmp(:,:) = Tracer%field(:,:,k,tau)*mass_t(:,:) tbar = tbar + mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = Tracer%field(:,:,k,tau)*Tracer%field(:,:,k,tau)*mass_t(:,:) tvar = tvar + mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) tmp(:,:) = abs(Tracer%field(:,:,k,taup1)-Tracer%field(:,:,k,taum1))*dtimer tchg = tchg + mpp_global_sum(Dom%domain2d,tmp(:,:), global_sum_flag) enddo 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 call mpp_max(num) 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 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 (material_flux(isd:ied,jsd:jed)) allocate (eta_smooth(isd:ied,jsd:jed)) allocate (pbot_smooth(isd:ied,jsd:jed)) material_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 material_flux(i,j) = material_flux(i,j) + dteta*Grd%tmask(i,j,1)*Grd%dat(i,j) & *( pme(i,j) + runoff(i,j) + calving(i,j) & + Ext_mode%source(i,j) + Ext_mode%eta_smooth(i,j) + Ext_mode%pbot_smooth(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) material_flux(:,:) = material_flux(:,:)*Grd%obc_tmask(:,:) if(have_obc) eta_smooth(:,:) = eta_smooth(:,:)*Grd%obc_tmask(:,:) if(have_obc) pbot_smooth(:,:) = pbot_smooth(:,:)*Grd%obc_tmask(:,:) mass_input = mpp_global_sum(Dom%domain2d,material_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) 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)') & ' Mass input via boundary fluxes over time interval= ',mass_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,'(a,es24.17,a)') & ' Change in ocean mass over time interval = ',mass_chg,' kg.' write (stdoutunit,'(a,es10.3,a)') & ' Error in mass content change = ',mass_error,' kg.' write (stdoutunit,'(a,es10.3,a)') & ' Error in rate of mass content change = ',mass_error_rate,' kg/s.' 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) :: runoff_mod real :: tracer_chg, tracer_error, tracer_error_rate real :: tracer_source_input, tracer_flux_input real :: tracer_eta_smooth_input, 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_flux(isd:ied,jsd:jed,num_prog_tracers)) 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_flux(:,:,:) = 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 ! include runoff and calving and pme contributions here ! (must remove from tracer_source below) do j=jsc,jec do i=isc,iec tracer_flux(i,j,n) = tracer_flux(i,j,n) + dtts*T_prog(n)%conversion*Grd%tmask(i,j,1)*Grd%dat(i,j) & *(-T_prog(n)%btf(i,j) + T_prog(n)%stf(i,j) + & T_prog(n)%tpme(i,j)*pme(i,j) + T_prog(n)%trunoff(i,j)*runoff_mod(i,j) + T_prog(n)%tcalving(i,j)*calving(i,j)) enddo enddo ! global integral of th_tendency leaves just the sources ! (and other stuff to be removed below) do k=1,nk do j=jsc,jec do i=isc,iec tracer_source(i,j,n) = tracer_source(i,j,n) + & dtts*T_prog(n)%conversion*Grd%tmask(i,j,k)*Grd%dat(i,j)*T_prog(n)%th_tendency(i,j,k) & +dtts*T_prog(n)%conversion*Grd%tmask(i,j,k)*Grd%dat(i,j)*Thickness%rho_dzt(i,j,k,taup1) & *T_prog(n)%source(i,j,k) enddo enddo enddo do j=jsc,jec do i=isc,iec tracer_eta_smooth(i,j,n) = tracer_eta_smooth(i,j,n) + & dtts*T_prog(n)%conversion*Grd%tmask(i,j,1)*Grd%dat(i,j)*T_prog(n)%eta_smooth(i,j) enddo enddo do j=jsc,jec do i=isc,iec tracer_pbot_smooth(i,j,n) = tracer_pbot_smooth(i,j,n) + & dtts*T_prog(n)%conversion*Grd%tmask(i,j,1)*Grd%dat(i,j)*T_prog(n)%pbot_smooth(i,j) enddo enddo ! remove runoff, calving, pme, and smooth from tracer source ! (they are added to T_prog(n)%th_tendency inside ocean_rivermix_mod and ocean_tracers) k=1 do j=jsc,jec do i=isc,iec tracer_source(i,j,n) = tracer_source(i,j,n) & -dtts*T_prog(n)%conversion*Grd%tmask(i,j,k)*Grd%dat(i,j) & *( T_prog(n)%trunoff(i,j)*runoff_mod(i,j) & +T_prog(n)%tcalving(i,j)*calving(i,j) & +T_prog(n)%tpme(i,j)*pme(i,j)) & -tracer_eta_smooth(i,j,n)-tracer_pbot_smooth(i,j,n) enddo enddo if(n==index_temp .and. index_frazil > 0) then do j=jsc,jec do i=isc,iec tracer_flux(i,j,n) = tracer_flux(i,j,n) & + T_diag(index_frazil)%field(i,j,1)*Grd%dat(i,j)*Grd%tmask(i,j,1) enddo enddo endif if (have_obc) then call ocean_obc_tracer_flux(Time, T_prog(n), tmp, n, send_out = .false.) tracer_flux(:,:,n) = (tracer_flux(:,:,n) + dtts*T_prog(n)%conversion*tmp(:,:))*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(:,:) 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_flux_input = mpp_global_sum(Dom%domain2d,tracer_flux(:,:,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_chg = tracer_final(n)-tracer_start(n) tracer_error = tracer_chg-tracer_flux_input-tracer_source_input & -tracer_eta_smooth_input-tracer_pbot_smooth_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)') ' Heat input by eta_t smoother over time interval = ',tracer_eta_smooth_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by pbot_t smoother over time interval = ',tracer_pbot_smooth_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input by sources over time interval = ',tracer_source_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Heat input through boundaries over time interval = ',tracer_flux_input,' J.' write (stdoutunit,'(a,es24.17,a)') ' Change in ocean heat content over time interval = ',tracer_chg,' J.' write (stdoutunit,'(a,es10.3,a)') ' Error in heat content change =',tracer_error,' J.' write (stdoutunit,'(a,es10.3,a)') ' Error in rate of heat content change =',tracer_error_rate,' W/m^2.' 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)// & ' input by eta_t smoother over time interval = ',tracer_eta_smooth_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by pbot_t smoother over time interval = ',tracer_pbot_smooth_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by sources over time interval = ',tracer_source_input,' kg*yr.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input through boundaries over time interval = ',tracer_flux_input,' kg*yr.' write (stdoutunit,'(a,es24.17,a)') ' Change in ocean '//trim(T_prog(n)%name)// & ' over time interval = ',tracer_chg,' kg*yr.' write (stdoutunit,'(a,es10.3,a)') ' Error in '//trim(T_prog(n)%name)//' change =', & tracer_error,' kg*yr.' write (stdoutunit,'(a,es10.3,a)') ' Error in rate of '//trim(T_prog(n)%name)//' change =', & sec_in_yr*tracer_error_rate,' kg/(m^2)' 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)// & ' input by eta_t smoother over time interval = ',tracer_eta_smooth_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by pbot_t smoother over time interval = ',tracer_pbot_smooth_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input by sources over time interval = ',tracer_source_input,' kg.' write (stdoutunit,'(1x,a,es24.17,a)') trim(T_prog(n)%name)// & ' input through boundaries over time interval = ',tracer_flux_input,' kg.' write (stdoutunit,'(a,es24.17,a)') ' Change in ocean '//trim(T_prog(n)%name)// & ' mass over time interval = ',tracer_chg,' kg.' write (stdoutunit,'(a,es10.3,a)') ' Error in '//trim(T_prog(n)%name)//' mass change =', & tracer_error,' kg' write (stdoutunit,'(a,es10.3,a)') ' Error in rate of '//trim(T_prog(n)%name)//' mass change =', & tracer_error_rate,' kg/(m^2 sec)' 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@noaa.gov ! ! ! 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 ! decreasing with depth. ! ! Author: Harper.Simmons@noaa.gov ! Zhi.Liang@noaa.gov ! ! 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 with depth. ! ! Based on "diagnose_depth_of_potrho" by Harper.Simmons@noaa.gov ! ! Author: Stephen.Griffies@noaa.gov ! Zhi.Liang@noaa.gov ! ! 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@noaa.gov ! ! 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 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_rho=1,potrho_nk do k=1,nk-1 do j=jsc,jec 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@noaa.gov ! 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 (Time, Thickness, Dens, potrho_mix_depth, potrho_mix_base) 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(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 (Time, 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 (Time, Thickness, Dens, & potrho_mix_depth = potrho_mix_depth) elseif (id_potrho_mix_base > 0) then call calc_potrho_mixed_layer (Time, Thickness, Dens, & potrho_mix_base = potrho_mix_base) endif if (id_potrho_mix_depth > 0) then used = send_data (id_potrho_mix_depth, potrho_mix_depth(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif if (id_potrho_mix_base > 0) then used = send_data (id_potrho_mix_base, potrho_mix_base(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif 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 :: i,j,k real :: tvolume 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) 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 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" end module ocean_tracer_diag_mod